Hello everyone,

I am estimating a DiD (difference in difference) with Least-squares dummy-variables (LSDV).

My data set contains 8,232 students (i) from 60 schools (j) and 460 classes (c) with a panel data format with T=5 (wave). For each student, I have the test scores (profic_mat) and a list of observed variables ($controlvar) over the time period.
During the time series (2003-2008), a policy change is implemented in state schools in year 2007. Then, students from state schools are my treatment group and students from municipal schools are the control group. My DiD is 1 if student is enrolled in state schools (treated) in post-treatment period (time).

Since the educational achievements follow a hierarchical structure, I assume individual, time, school and class fixed effects.

Fitting two-way fixed-effects models is relatively simple. I estimate a model with individual and time fixed effects including the time effects as dummies and eliminating the individual effects by the withhin transformation.

Code:
xtreg profic_mat DiD time treated $controlvar wave_*, fe i(IDstudent) nonest cluster(IDschool)
By defining a spell, I can fit a three-way fixed-effects. Then I treat as one spell all unique combinations of i and j.
Code:
egen spell=group(IDstudent IDschool)
xtreg profic_mat DiD time treated $controlvar wave_*, fe i(spell) cluster(IDschool)
My question now is how can I fit a model with a four-way error-components?
I have estimated the two following models, but no one produced reliable values.

Code:
egen spellThree=group(IDstudent IDschool IDclass)
xtreg profic_mat DiD time treated $controlvar wave_*, fe i(spellThree) cluster(IDschool)
                                                           /*(AND)*/
xtreg profic_mat DiD time treated $controlvar wave_* IDclass_*, fe i(IDstudent) nonest cluster(IDschool)
Any advice would be highly appreciated!