load x[200,5]=grunfeld.m;
y=x[.,3];
xs=x[.,4:5];

n=10;
t=20;
nt=n*t;

jnt=ones(200,1);
jt=ones(t,1);
jn=ones(n,1);
xp=(jnt~xs);
bp=invpd(xp'xp)*xp'*y;
"pooled :" bp;

D=eye(n) .*. jt; 
pd=D*invpd(d'd)*d';
md=eye(nt)-pd;

y1=md*y;
x1=md*xs;

bfe=invpd(x1'x1)*x1'*y1;
"fixed " bfe;

x3=pd*xp;
bbe2=invpd(x3'x3)*x3'y2;
"between " bbe2;

x4=(eye(n).*.jt)~xs;
bfe1=invpd(x4'*x4)*x4'y;
"Fixed 2 :" bfe1;
ac = meanc(bfe1[1:n]);
"average constant : " ac;