/* STEINUR.PRG Calculations reported in "Asymptotic moments of autoregressive estimates with a near unit root" by Bruce E. Hansen Department of Economics, University of Wisconsin This is a Gauss program */ declare p,c; library pgraph; graphset; pqgwin many; _pdate=""; _plwidth=10; let _pmcolor = 0 0 0 0 0 0 0 0 15; fonts("complex simgrma"); _pcolor=1; _ptitlht=.15; inc1=1; inc2=20; cmax1=20; cmax2=400; @ Grid on local-to-unity parameter c @ cn1=cmax1/inc1+1; cn2=(cmax2-cmax1)/inc2; cn=cn1+cn2; cs=seqa(0,-inc1,cn1)|seqa(-cmax1-inc2,-inc2,cn2); @ Calculating Moments @ gr=100; @ NUmber of gridpoints on [0,1] @ _intord=40; @ Order of Gauss-Legendre quadrature @ m=zeros(cn,4); xx=(seqa(1/gr,1/gr,gr)~seqa(0,1/gr,gr))'; for pp (1,4,1); p=pp; for ci (1,cn,1); c=cs[ci]; intq=intquad1(&g,xx); m[ci,pp]=sumc(intq); endfor; endfor; @ mu_r(c) @ m1=m[.,1]-cs; m2=m[.,2]-m[.,1].*cs*2+(cs.^2); m3=m[.,3]-m[.,2].*cs*3+m[.,1].*(cs.^2)*3-(cs.^3); m4=m[.,4]-m[.,3].*cs*4+m[.,2].*(cs.^2)*6-m[.,1].*(cs.^3)*4+(cs.^4); @ Moments @ mse=m2; var=m2-(m1.^2); skew=(m3-m2.*m1*3+(m1.^3)*2)./(var.^(3/2)); kur=(m4-(m3.*m1)*4+(m2.*m1.*m1)*6-(m1.^4)*3)./(var.^2); c1=cs+(cs.==0); gc=-(1+(1-exp(c1*2))./(c1*2))./(c1*2); gc=gc.*(cs .ne 0) +(cs.==0)/2; stmean=sqrt(gc).*m1; stvar=gc.*var; stmse=gc.*mse; @ Simulation @ rep=1000000; let ns = 500 1000 10000 100000; nn=rows(ns); means=zeros(cn,nn); mses=zeros(cn,nn); vars=zeros(cn,nn); skews=zeros(cn,nn); kurs=zeros(cn,nn); mse2=zeros(cn,nn); for ni (1,nn,1); n=ns[ni]; a=1+cs'/n; mu1=zeros(cn,1); mu2=zeros(cn,1); mu3=zeros(cn,1); mu4=zeros(cn,1); mu2s=zeros(cn,1); for i (1,rep,1); e=rndn(n,cn); y=recserar(e,e[1,.],a); vc=meanc(y.^2)/n; uc=((y[n,.].^2)'/n-1)/2; uv=uc./vc; w=1-vc./(uc.^2); w=w.*(w.>0); uvs=uv.*w.*(w.>0); uvs=uvs-cs; uv=uv-cs; mu1=mu1+uv; mu2=mu2+(uv.^2); mu3=mu3+(uv.^3); mu4=mu4+(uv.^4); mu2s=mu2s+(uvs.^2); endfor; mu1=mu1./rep; mu2=mu2./rep; mu3=mu3./rep; mu4=mu4./rep; mu2s=mu2s./rep; means[.,ni]=mu1; mses[.,ni]=mu2; mse2[.,ni]=mu2s; vars[.,ni]=mu2-(mu1.^2); skews[.,ni]=(mu3-mu2.*mu1*3+(mu1.^3)*2)./((mu2-(mu1.^2)).^(3/2)); kurs[.,ni]=(mu4-(mu3.*mu1)*4+(mu2.*mu1.*mu1)*6-(mu1.^4)*3)./((mu2-(mu1.^2)).^2); endfor; stmses=gc.*mses; stmse2=gc.*mse2; output file=int.out reset; "Replications = " rep; "Sample size = " ns[nn]; "Max R(Stein)-R(OLS) = " maxc(mse2[.,nn]-mses[.,nn]); "Max R(OLS)-R(Stein) = " maxc(mses[.,nn]-mse2[.,nn]); "Max SR(Stein)-R(OLS) = " maxc(stmse2[.,nn]-stmses[.,nn]); "Max SR(OLS)-R(Stein) = " maxc(stmses[.,nn]-stmse2[.,nn]); output off; output file=steinur.out reset; cs~m1~var~skew~kur; output off; _plegstr= "" \$+ "Exact"; for ni (1,nn,1); _plegstr=_plegstr \$+ "\000n=" \$+ ftocv(ns[ni],1,0); endfor; c1=cs[1:cn1]; xtics(-cmax1,0,2,2); @title("Mean");@ ytics(-2.2,-1.2,.2,10); let _plegctl= 1 5 -5 -1.5; xy(c1,m1[1:cn1,.]~means[1:cn1,.]); @title("Variance");@ ytics(0,50,5,10); let _plegctl= 1 5 -5 37; xy(c1,var[1:cn1,.]~vars[1:cn1,.]); @title("Skewness");@ ytics(-2.5,0,.5,10); let _plegctl= 1 5 -5 -0.7; xy(c1,skew[1:cn1,.]~skews[1:cn1,.]); @title("Kurtosis");@ ytics(4,14,2,10); let _plegctl= 1 5 -20 11.6; xy(c1,kur[1:cn1,.]~kurs[1:cn1,.]); _plegstr= "Stein\000OLS"; _plegctl=1; rmse=mse2[.,nn]./mses[.,nn]; r1=ones(cn,1); @title("MSE");@ ytics(0,1.2,.2,10); xy(c1,rmse[1:cn1]~r1[1:cn1]); _plegstr= "" \$+ "Exact"; for ni (1,nn,1); _plegstr=_plegstr \$+ "\000n=" \$+ ftocv(ns[ni],1,0); endfor; xtics(-cmax2,0,100,5); @title("Mean");@ ytics(-20,140,40,10); let _plegctl= 1 5 -100 100; xy(cs,m1~means); @title("Variance");@ ytics(0,800,100,10); let _plegctl= 1 5 -100 620; xy(cs,var~vars); @title("Skewness");@ ytics(-3,0,1,10); let _plegctl= 1 5 -400 -3; xy(cs,skew~skews); @title("Kurtosis");@ ytics(2,12,1,10); let _plegctl= 1 5 -400 9.7; xy(cs,kur~kurs); _plegstr= "Stein\000OLS"; _plegctl=1; ytics(0,1.2,.2,10); xy(cs,rmse~r1); /***********************************************/ /* Takes p and c as given (globals) */ proc g(x); local x1,lx,lx0,psi,f,a,aa,cst; x1=x.*(x.<1); lx=x./(1-x1) - c; psi=gpsi(lx); f=lx.*(((lx.^2)-(c^2)).^(p-1))./((1-x1).^2).*exp(-(x./(1-x1))/2)./sqrt(1+exp(-lx*2)); f=f./sqrt(1-psi*c); a=(-1)^p; aa=a; for j (1,p,1); a=-(a.*psi./(1-psi*c))*(j*2-1); aa=aa+a*numcombinations(p,j); endfor; f=f.*aa; f=f.*(x.<1); cst=((p-1)!)*(2^(2*p-3/2)); f=f./cst; retp(f); endp; @ Robust version of psi(u)=tanh(u)/u @ proc gpsi(x); local t1,t2,x1,tx,y,ax; ax=abs(x); tx=tanh(ax.*(ax .<= 700))+(ax .> 700); t1=(ax .< 1e-4); t2=1-t1; x1=ax.*t2+t1; y=t1+(tx./x1).*t2; retp(y); endp;