*
* Example 16.6 from page 449
*
open data tablef9-2[1].txt
data(format=prn,org=columns) 1 25 valueadd capital labor nfirm
*
set logy = log(valueadd)
set logk = log(capital)
set logl = log(labor)
*
linreg logy / resols
# constant logk logl
linreg(smpl=t<>4.and.t<>10) logy
# constant logk logl
rreg logy / reslad
# constant logk logl
compute betalad=%beta
*
prj fitted
compute [vector] bboot=%zeros(%nreg,1)
compute [symm]   bsq  =%zeros(%nreg,%nreg)
do i=1,5000
   boot entries
   set yx = logy(entries(t))
   set kx = logk(entries(t))
   set lx = logl(entries(t))
   rreg(noprint) yx
   # constant kx lx
   compute bsq  =bsq  +%outerxx(%beta-betalad)
end do i
disp "Bootstrapped standard errors" %sqrt(%xdiag(bsq*(1.0/5000)))
*
compute [vector] bboot=%zeros(%nreg,1)
compute [symm]   bsq  =%zeros(%nreg,%nreg)
do i=1,5000
   boot entries
   set yx = fitted+reslad(entries(t))
   rreg(noprint) yx
   # constant logk logl
   compute bsq  =bsq  +%outerxx(%beta-betalad)
end do i
disp "Bootstrapped standard errors" %sqrt(%xdiag(bsq*(1.0/5000)))
*
* Bayesian estimates
* Use adjusted OLS standard errors with the LAD estimator as an importance
* function.
*
linreg(noprint) logy
# constant logk logl
compute xxs=%decomp(%pi*.5*%seesq*%xx)
rreg(noprint) logy / reslad
# constant logk logl
set absu = abs(reslad)
compute sumabsu=%sum(absu)
compute lambda=(1.0/%nobs)*sumabsu
dec vector beta u(%nreg)
compute betalad=%beta
compute [vector] bbayes=%zeros(%nreg,1)
compute [symm]   bayesq=%zeros(%nreg,%nreg)
compute wsum=0.0
do i=1,500
   compute u=%ran(1.0)
   compute [vector] beta=betalad+xxs*u
   set absu = abs(logy-beta(1)-beta(2)*logk-beta(3)*logl)
   compute sumabsu=%sum(absu)
   compute w=exp(-sumabsu/lambda+%nobs+.5*%dot(u,u))
   compute bbayes=bbayes+w*beta
   compute bayesq=bayesq+w*%outerxx(beta)
   compute wsum=wsum+w
end do i
*
disp "Bayesian mean" bbayes*(1.0/wsum)
ewise bayesq(i,j)=bayesq(i,j)/wsum-bbayes(i)*bbayes(j)/(wsum**2)
disp "Bayesian std errors" %sqrt(%xdiag(bayesq))