alpha=-1.342511; beta=.9914399;
load lphdata
z=[ones(111,1),hsdata(2:112,1)./hsdata(1:111,1),(hsdata(2:112,2)+ ..
   hsdata(2:112,3))./hsdata(1:111,2)];

h=beta*z(:,3).*z(:,2).^alpha-1;
dhdb=[(h+1)/beta,log(z(:,2)).*z(:,3).*z(:,2).^alpha];
z=z(1:110,:);
tem=z'*h(2:111)*h(2:111)'*z
f=h(2:111)'*z*inv(z'*z)*z'*h(2:111);
cov=dhdb(2:111,:)'*z*inv(z'*z)*z'*dhdb(2:111,:)
