/* 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;