/*================================================
Panel Analysis : Fixed Effect Model

Yit = alpha_i + beta*X_it  +  E_it

E_it ~ i.i.d N(0, sig2)
i = 1, . . , n
t = 1, 2, . . ,T
================================================*/

new;
cls;

/*==================== DGP =======================*/
t = 30;    @  size of times series @
n = 3;     @ # of Cross section @
k = 1;     @ # of regressors @

/*=== True parameters===*/
mu1 = 2;
mu2 = 0;
mu3 = -2;
tru_b = 1;
sig2 = 2;

tru_para = tru_b|sig2;

x1 = rndu(T,1)*5+mu1;
x2 = rndu(T,1)*5+mu2;
x3 = rndu(T,1)*5+mu3;
x = x1|x2|x3;

y1 = mu1 + x1*tru_b+ rndn(T,1)*sqrt(sig2);
y2 = mu2 + x2*tru_b+ rndn(T,1)*sqrt(sig2);
y3 = mu3 + x3*tru_b+ rndn(T,1)*sqrt(sig2);

y = y1|y2|y3;

/*========== Estimation 1: Fixed Effect Estimator ============*/
@ Demean for individual 1 @
y1_dm = y1 - meanc(y1);
x1_dm = x1 - meanc(x1);

@ Demean for individual 2 @
y2_dm = y2 - meanc(y2);
x2_dm = x2 - meanc(x2);

@ Demean for individual 3 @
y3_dm = y3 - meanc(y3);
x3_dm = x3 - meanc(x3);

y_dm = y1_dm|y2_dm|y3_dm;
x_dm = x1_dm|x2_dm|x3_dm;

b_hat = inv(x_dm'x_dm)*x_dm'y_dm; @ Fixed Effect Estimator @
ehat = y_dm-x_dm*b_hat;
sig2hat = ehat'ehat/(n*T-n-k);  @ n=3, k=1 @
cov = sig2hat*inv(x_dm'x_dm);

"=============<Results>=====================";
"True value                     ";;tru_para';
"Fixed Effect Estimates  ";;b_hat'~sig2hat;
"S.D of Estimates            ";;sqrt(diag(cov))';
"==========================================";

/*========== Estimation 2: The Pooled Estimator ============*/
Ymat = y;
Xmat = ones(n*T,1)~X;
B_pool = inv(Xmat'Xmat)*Xmat'Ymat;
"==========================================";
"Pooled Estimates          ";;b_pool[2];
"==========================================";


end;