load lphdata
x=[hsdata(2:112,1)./hsdata(1:111,1),(hsdata(2:112,2)+ ..
   hsdata(2:112,3))./hsdata(1:111,2)];
z=[ones(105,1),x(1:105,:),x(2:106,:),x(3:107,:),x(4:108,:),x(5:109,:),x(6:110,:)];
x=x(7:111,:);
q=13; r=13; l=2; T=105;

alpha=-3.436630; beta=1.006194;
H1=beta*x(:,2).*x(:,1).^alpha-1;
dH1=[(H1+1)/beta,log(x(:,1)).*x(:,2).*x(:,1).^alpha];
F1=(ones(q,1)*H1').*z';
g1=sum(F1')'/T;
W1=inv(z'*z/T);
J1=g1'*W1*g1
dF1=[(ones(q,1)*dH1(:,1)').*z';(ones(q,1)*dH1(:,2)').*z'];
D1=zeros(r,l);
tem=sum(dF1')'/T;
D1(:)=tem;
vcov1=inv(D1'*W1*D1);
se1=sqrt(diag(vcov1))

alpha=-3.432791; beta=1.005945;
H2=beta*x(:,2).*x(:,1).^alpha-1;
dH2=[(H2+1)/beta,log(x(:,1)).*x(:,2).*x(:,1).^alpha];
F2=(ones(q,1)*H2').*z';
g2=sum(F2')'/T;
R2=F2*F2'/T;
W2=inv(R2);
J2=g2'*W2*g2
dF2=[(ones(q,1)*dH2(:,1)').*z';(ones(q,1)*dH2(:,2)').*z'];
D2=zeros(r,l);
tem=sum(dF2')'/T;
D2(:)=tem;
vcov2=inv(D2'*W2*D2);
se2=sqrt(diag(vcov2))
