This code constructs the design matrix for a panel with both time and individual fixed effects and then estimates the model.


\[Beta] = 3;
NumIndividuals = 323;
NumPeriods = 4;
NumObs = NumIndividuals * NumPeriods;
x = Table[Random[NormalDistribution[0, 1]], {NumObs}];
y = \[Beta]*x + Table[Random[NormalDistribution[0, 1]], {NumObs}];
XindivFull =
KroneckerProduct[IdentityMatrix[NumPeriods],
Table[1, {NumIndividuals}]] // Transpose;
XtimeFull =
KroneckerProduct[Table[1, {NumIndividuals}],
IdentityMatrix[NumPeriods]];
Xcombo = Join[XindivFull[[All, {2, NumPeriods}]],
XtimeFull[[All, {2, NumPeriods}]], 2];
\[Iota] = Table[1, {NumObs}];
XwIntercept = MapThread[Append, {Xcombo, \[Iota]}];
Xfull = MapThread[Append, {XwIntercept, x}];
\[Beta]hat = Inverse[(Transpose[Xfull].Xfull)].Transpose[Xfull].y