-------------------------------------------------------------------------------
name: <unnamed>
log: /Users/sam/statadev/mcmclinear_all/mcmclinear_dev/mcmclinear.log
log type: text
opened on: 5 Jan 2012, 10:11:25
.
. version 12.0
.
. mata: mata mlib create lmcmclinear, replace
(file lmcmclinear.mlib created)
.
. foreach f in mcmclinear_reg mcmclinear_mixed {
2. do `f'.mata
3. mata: mata mlib add lmcmclinear `f'()
4. }
. version 12.0
. mata:
------------------------------------------------- mata (type end to exit) -----
: mata set matastrict on
: mata set matafavor speed
:
: /* version 1.0, 4 jan 2012 */
: /* sam schulhofer-wohl, federal reserve bank of minneapolis */
:
: /* gibbs sampling for linear regression */
: void mcmclinear_reg(string scalar yy, string scalar xx,
> real scalar iters, real scalar seed, real scalar d0,
> real scalar useconstant, string scalar ww, real scalar useweight,
> real scalar weight_is_aweight, real scalar showiter, real scalar usex)
> {
>
> /* declarations */
> real matrix y, x, w, wsqrt, xPx, xPy, betahat, xPx_ci, myP2x, beta, beta_ou
> t, sigma2_out
> real scalar Kx, lx, df, yPy, sigma, invsig2, iter, i, sigma2
>
> /* set seed */
> rseed(seed)
>
> /* find data */
> y=st_data(.,yy)
> if(usex) x=st_data(.,xx)
> if(useweight) w=st_data(.,ww)
>
> /* transform data if we have aweights */
> if(weight_is_aweight & useweight) {
> wsqrt=sqrt(w)
> y=y:*wsqrt
> if(useconstant) {
> if(usex) x=(x,J(rows(x),1,1)):*wsqrt
> else {
> x=wsqrt
> usex=1
> }
> useconstant=0
> }
> else x=x:*wsqrt /* always have either useconstant or usex */
> useweight=0
> }
>
> /* useful scalars */
> if(usex) Kx=cols(x)
> else Kx=0
> if(useconstant) lx=Kx+1
> else lx=Kx
>
> /* precalculate matrices used repeatedly */
> if(useweight) {
> if(usex) {
> xPx=quadcross(x,useconstant,w,x,useconstant)
> xPy=quadcross(x,useconstant,w,y,0)
> }
> else {
> xPx=quadsum(w)
> xPy=quadcross(w,y)
> }
> yPy=quadcross(y,w,y)
> }
> else {
> if(usex) {
> xPx=quadcross(x,useconstant,x,useconstant)
> xPy=quadcross(x,useconstant,y,0)
> }
> else {
> xPx=rows(y)
> xPy=quadsum(y)
> }
> yPy=quadcross(y,y)
> }
> xPx_ci=cholesky(luinv(xPx))
> myP2x=-2*xPy'
> betahat=lusolve(xPx,xPy)
> if(useweight) df=(quadsum(w)+1)/2
> else df=(rows(y)+1)/2
>
> /* set up initial guess for sigma2 */
> /* changing seed will change this initial guess, so you can get different c
> hains */
> /* this is the only relevant initial guess */
> invsig2=rgamma(1,1,df,2/(yPy+(myP2x+betahat'*xPx)*betahat+d0))
> sigma2=1/invsig2
> sigma=sqrt(sigma2)
>
> /* matrices to store results */
> beta_out=(betahat'\J(iters,lx,0))
> sigma2_out=(sigma2\J(iters,1,0))
>
> /* run the gibbs sampler */
> for(iter=1;iter<=iters;iter++) {
>
> if(showiter) {
> display("iteration "+strofreal(iter))
> displayflush()
> }
>
> /* 1. draw beta */
> beta=betahat+sigma*xPx_ci*rnormal(lx,1,0,1)
> beta_out[iter+1,.]=beta'
>
> /* 2. draw sigma2 */
> invsig2=rgamma(1,1,df,2/(yPy+(myP2x+betahat'*xPx)*betahat+d0))
> sigma2=1/invsig2
> sigma=sqrt(sigma2)
> sigma2_out[iter+1,1]=sigma2
>
> }
>
> /* return results */
> stata("drop _all")
> st_addobs(iters+1)
> for(i=1;i<=lx;i++) {
> (void) st_addvar("double","beta"+strofreal(i))
> }
> st_store(.,range(1,lx,1)',beta_out)
> (void) st_addvar("double","sigma2")
> st_store(.,lx+1,sigma2_out)
>
> }
:
: end
-------------------------------------------------------------------------------
.
end of do-file
(1 function added)
. version 12.0
. mata:
------------------------------------------------- mata (type end to exit) -----
: mata set matastrict on
: mata set matafavor speed
:
: /* gibbs sampling for linear mixed model */
: void mcmclinear_mixed(string scalar yy, string scalar xx, string scalar zz,
> string scalar gg, real scalar iters, real scalar seed, real scalar d0, stri
> ng scalar ddelta0,
> real scalar useconstant_fe, real scalar useconstant_re,
> string scalar ww_fe, real scalar useweight_fe,
> real scalar feweight_is_aweight,
> real scalar showiter, real scalar N_g, string scalar sigrownames, string sc
> alar sigcolnames,
> real scalar usex, real scalar usez)
> {
>
> /* declarations */
> real matrix y, x, z, g, w_fe, wsqrt_fe, beta, beta_out, sigma2_out, Sigma_o
> ut, theta_g_out, ginfo,
> xPx_g, zPz_g, xPy_g, xPz_g, zPxy_g, zPy_g, zPx_g, m2yPx_g, m2yPz_g, xP2z_
> g,
> work_x, work_y, work_z, work_w, invsigma2_zPz_g, invsigma2_zPxy_g, theta_
> g, Sigma, work, D, Dy, E, Sigmainv, delta0
> real scalar Kx, lx, df, sigma2, invsig2, iter, i, gi, Kz, lz, chi_sig2eps_s
> tart, chi_sig2eps, usecfetemp, usecretemp
>
> /* set seed */
> rseed(seed)
>
> /* find data */
> y=st_data(.,yy)
> if(usex) x=st_data(.,xx)
> if(usez) z=st_data(.,zz)
> g=st_data(.,gg)
> if(useweight_fe) w_fe=st_data(.,ww_fe)
> delta0=st_matrix(ddelta0)
>
> /* set up panel identifiers */
> ginfo=panelsetup(g,1)
>
> /* transform data if we have aweights */
> if(feweight_is_aweight & useweight_fe) {
> wsqrt_fe=sqrt(w_fe)
> y=y:*wsqrt_fe
> if(useconstant_fe) {
> if(usex) x=(x,J(rows(x),1,1)):*wsqrt_fe
> else {
> x=wsqrt_fe
> usex=1
> }
> useconstant_fe=0
> }
> else x=x:*wsqrt_fe /* always have either useconstant_fe or usex */
> useweight_fe=0
> }
>
> /* useful scalars */
> if(usex) Kx=cols(x)
> else Kx=0
> if(useconstant_fe) lx=Kx+1
> else lx=Kx
> if(usez) Kz=cols(z)
> else Kz=0
> if(useconstant_re) lz=Kz+1
> else lz=Kz
>
> /* precalculate matrices used repeatedly */
> xPx_g=J(lx,lx*N_g,.)
> zPz_g=J(lz,lz*N_g,.)
> xPy_g=J(lx,N_g,.)
> xPz_g=J(lx,lz*N_g,.)
> zPy_g=J(lz,N_g,.)
> zPx_g=J(lz,lx*N_g,.)
> zPxy_g=J(lz,(lx+1)*N_g,.)
> m2yPx_g=J(N_g,lx,.)
> m2yPz_g=J(N_g,lz,.)
> xP2z_g=J(lx,N_g*lz,.)
> if(usex) usecfetemp=useconstant_fe
> else usecfetemp=0
> if(usez) usecretemp=useconstant_re
> else usecretemp=0
> for(gi=1;gi<=N_g;gi++) {
> work_y=y[|ginfo[gi,1],1\ginfo[gi,2],1|]
> if(usex) work_x=x[|ginfo[gi,1],1\ginfo[gi,2],.|]
> else work_x=J(rows(work_y),1,1)
> if(usez) work_z=z[|ginfo[gi,1],1\ginfo[gi,2],.|]
> else work_z=J(rows(work_y),1,1)
> if(useweight_fe) {
> work_w=w_fe[|ginfo[gi,1],1\ginfo[gi,2],1|]
> xPx_g[|1,(gi-1)*lx+1\.,gi*lx|]=quadcross(work_x, usecfetemp,work_w,work
> _x, usecfetemp)
> xPy_g[.,gi]=quadcross(work_x, usecfetemp,work_w,work_y,0)
> xPz_g[|1,(gi-1)*lz+1\.,gi*lz|]=quadcross(work_x, usecfetemp,work_w,work
> _z, usecretemp)
> zPy_g[.,gi]=quadcross(work_z, usecretemp,work_w,work_y,0)
> zPz_g[|1,(gi-1)*lz+1\.,gi*lz|]=quadcross(work_z, usecretemp,work_w,work
> _z, usecretemp)
> }
> else {
> xPx_g[|1,(gi-1)*lx+1\.,gi*lx|]=quadcross(work_x, usecfetemp,work_x, use
> cfetemp)
> xPy_g[.,gi]=quadcross(work_x, usecfetemp,work_y,0)
> xPz_g[|1,(gi-1)*lz+1\.,gi*lz|]=quadcross(work_x, usecfetemp,work_z, use
> cretemp)
> zPy_g[.,gi]=quadcross(work_z, usecretemp,work_y,0)
> zPz_g[|1,(gi-1)*lz+1\.,gi*lz|]=quadcross(work_z, usecretemp,work_z, use
> cretemp)
> }
> zPx_g[|1,(gi-1)*lx+1\.,gi*lx|]=xPz_g[|1,(gi-1)*lz+1\.,gi*lz|]'
> zPxy_g[|1,(gi-1)*(lx+1)+1\.,gi*(lx+1)|]=(zPx_g[|1,(gi-1)*lx+1\.,gi*lx|],z
> Py_g[.,gi])
> }
> m2yPx_g=-2:*xPy_g'
> m2yPz_g=-2:*zPy_g'
> xP2z_g=2:*xPz_g
> if(useweight_fe) chi_sig2eps_start=d0+quadcross(y,w_fe,y)
> else chi_sig2eps_start=d0+quadcross(y,y)
> if(useweight_fe) df=(quadsum(w_fe)+1)/2
> else df=(rows(y)+1)/2
>
> /* set up initial guesses for sigma2 and Sigma */
> /* strategy here is to guess beta and theta_g by running regressions, then
> draw random estimates of sigma2 and Sigma based on the guesses for beta and t
> heta_g */
> /* the guesses for sigma2 and Sigma are all that matter for initializing th
> e chain; thus, different seeds give you different chains */
> if(useweight_fe) {
> if(usex) beta=cholsolve(quadcross(x,useconstant_fe,w_fe,x,useconstant_fe)
> ,quadcross(x,useconstant_fe,w_fe,y,0))
> else beta=mean(y,w_fe)
> }
> else {
> if(usex) beta=cholsolve(quadcross(x,useconstant_fe,x,useconstant_fe),quad
> cross(x,useconstant_fe,y,0))
> else beta=mean(y)
> }
> theta_g=J(lz,N_g,0)
> for(gi=1;gi<=N_g;gi++) {
> E=cholinv(delta0)+(1/d0):*zPz_g[|1,(gi-1)*lz+1\.,gi*lz|]
> theta_g[.,gi]=cholsolve(E,(1/d0):*(zPxy_g[.,gi*(lx+1)]-zPxy_g[|1,(gi-1)*(
> lx+1)+1\.,gi*(lx+1)-1|]*beta))+cholesky(cholinv(E))*rnormal(lz,1,0,1)
> }
> work=(cholesky(luinv(delta0+quadcross(theta_g',theta_g')))*rnormal(lz,1+N_g
> ,0,1))'
> Sigmainv=quadcross(work,work)
> Sigma=luinv(Sigmainv)
> chi_sig2eps=chi_sig2eps_start
> work=J(1,lx,0)
> for(gi=1;gi<=N_g;gi++) {
> work=work+m2yPx_g[gi,.]+beta'*xPx_g[|1,(gi-1)*lx+1\.,gi*lx|]
> chi_sig2eps=chi_sig2eps+(m2yPz_g[gi,.]+beta'*xP2z_g[|1,(gi-1)*lz+1\.,gi*l
> z|]+theta_g[.,gi]'*zPz_g[|1,(gi-1)*lz+1\.,gi*lz|])*theta_g[.,gi]
> }
> invsig2=rgamma(1,1,df,2/(chi_sig2eps+work*beta))
> sigma2=1/invsig2
>
> /* matrices to store results */
> beta_out=(beta'\J(iters,lx,0))
> theta_g_out=(vec(theta_g)'\J(iters,lz*N_g,0))
> sigma2_out=(sigma2\J(iters,1,0))
> Sigma_out=(vech(Sigma)'\J(iters,lz*(lz+1)/2,0))
>
> /* run the gibbs sampler */
> for(iter=1;iter<=iters;iter++) {
>
> if(showiter) {
> display("iteration "+strofreal(iter))
> displayflush()
> }
>
> /* 1. matrices used repeatedly */
> invsigma2_zPz_g=invsig2:*zPz_g
> invsigma2_zPxy_g=invsig2:*zPxy_g
>
> /* 2. draw beta */
> D=J(lx,lx,0)
> Dy=J(lx,1,0)
> for(gi=1;gi<=N_g;gi++) {
> work=xPz_g[|1,(gi-1)*lz+1\.,gi*lz|]*cholsolve(Sigmainv+invsigma2_zPz_g[
> |1,(gi-1)*lz+1\.,gi*lz|],invsigma2_zPxy_g[|1,(gi-1)*(lx+1)+1\.,gi*(lx+1)|])
> D=D+xPx_g[|1,(gi-1)*lx+1\.,gi*lx|]-work[|1,1\.,lx|]
> Dy=Dy+xPy_g[.,gi]-work[.,lx+1]
> }
> beta=cholsolve(D,Dy)+cholesky(cholinv(D))*rnormal(lx,1,0,sqrt(sigma2))
> beta_out[iter+1,.]=beta'
>
> /* 3. draw theta_g */
> for(gi=1;gi<=N_g;gi++) {
> E=Sigmainv+invsigma2_zPz_g[|1,(gi-1)*lz+1\.,gi*lz|]
> theta_g[.,gi]=cholsolve(E,invsigma2_zPxy_g[.,gi*(lx+1)]-invsigma2_zPxy_
> g[|1,(gi-1)*(lx+1)+1\.,gi*(lx+1)-1|]*beta)+cholesky(cholinv(E))*rnormal(lz,1,
> 0,1)
> }
> theta_g_out[iter+1,.]=vec(theta_g)'
>
> /* 4. draw Sigma */
> work=(cholesky(luinv(delta0+quadcross(theta_g',theta_g')))*rnormal(lz,1+N
> _g,0,1))'
> Sigmainv=quadcross(work,work)
> Sigma=luinv(Sigmainv)
> Sigma_out[iter+1,.]=vech(Sigma)'
>
> /* 5. draw sigma2 */
> chi_sig2eps=chi_sig2eps_start
> work=J(1,lx,0)
> for(gi=1;gi<=N_g;gi++) {
> work=work+m2yPx_g[gi,.]+beta'*xPx_g[|1,(gi-1)*lx+1\.,gi*lx|]
> chi_sig2eps=chi_sig2eps+(m2yPz_g[gi,.]+beta'*xP2z_g[|1,(gi-1)*lz+1\.,gi
> *lz|]+theta_g[.,gi]'*zPz_g[|1,(gi-1)*lz+1\.,gi*lz|])*theta_g[.,gi]
> }
> invsig2=rgamma(1,1,df,2/(chi_sig2eps+work*beta))
> sigma2=1/invsig2
> sigma2_out[iter+1,1]=sigma2
>
> }
>
> /* return results */
> stata("drop _all")
> st_addobs(iters+1)
> for(i=1;i<=lx;i++) {
> (void) st_addvar("double","beta"+strofreal(i))
> }
> st_store(.,range(1,lx,1)',beta_out)
> (void) st_addvar("double","sigma2")
> st_store(.,lx+1,sigma2_out)
> for(i=1;i<=lz*N_g;i++) {
> (void) st_addvar("double","theta"+strofreal(i))
> }
> st_store(.,range(lx+2,lx+1+lz*N_g,1)',theta_g_out)
> for(i=1;i<=cols(Sigma_out);i++) {
> (void) st_addvar("double","Sigma"+strofreal(i))
> }
> st_store(.,range(lx+1+lz*N_g+1,lx+1+lz*N_g+cols(Sigma_out),1)',Sigma_out)
>
> st_matrix(sigrownames,vech(range(1,lz,1)#J(1,lz,1)))
> st_matrix(sigcolnames,vech(J(lz,1,1)#(range(1,lz,1)')))
>
> }
:
: end
-------------------------------------------------------------------------------
.
end of do-file
(1 function added)
.
. mata: mata mlib index
.mlib libraries to be searched are now
lmatabase;lmataado;lmataopt;lmatasem;lmcmclinear;l_cfrmt;lfreduse;lsmwoodbu
> ry;lxtabond2
.
.
. log close
name: <unnamed>
log: /Users/sam/statadev/mcmclinear_all/mcmclinear_dev/mcmclinear.log
log type: text
closed on: 5 Jan 2012, 10:11:25
-------------------------------------------------------------------------------