[学科前沿] Gauss 编程应用于学科前沿一：单位根 [推广有奖]

31楼

xuelida

发表于 2012-4-12 19:09:27 |只看作者 |坛友微信交流群

这是线性单位根检验的最后一例了，包括很多的unit root test的改进
Phillips-Perron Za
Modified Phillips-Perron MZa
Modified Phillips-Perron MZt
Modified Sargan Bhargava
ERS pt
Modified pt
DF-GLS
Said-Dickey-Fuller
MZa with GLS detrended data

运行程序：
new;cls;
bigt=100;
seed=100;
e=rndns(bigt,1,seed);
y=recserar(e,0,1.0); @生成自回归数据@
penalty=0;  /* 0 for maic 1 for bic*/ @0对应于maic，1对应于bic @
p=0;       /* order of deterministic terms */
kmin=0;
kmax=12*int(bigt/100)^(1/4);
{za,mza,msb,adf,pt,mpt,mzt,za1,mza1,mza2,dfols,krule,krule1,a1,a1ols}
=mic1(y,p,penalty,kmax,kmin);
print "k                         " krule;
print "alpha-hat                   " a1;
print "Tests Always Using GLS detrended data";
print "                            test    5% c.v.";
print "Phillips-Perron Za          " za';
print "Modified Phillips-Perron MZa " mza';
print "Modified Phillips-Perron MZt " mzt';
print "Modified Sargan Bhargava    " msb';
print "ERS pt                      " pt';
print "Modified pt                " mpt';
print "DF-GLS                      "  adf';
print;
print "Tests based on OLS detrending";
print "k       " krule1;
print "alpha-hat                   " a1ols;
print "Phillips-Perron Za          " za1';
print "Modified Phillips-Perron MZa " mza1';
print "Said-Dickey-Fuller          " dfols';
print "MZa with GLS detrended data; but s2ar is based on OLS detrending" mza2';

使用道具举报

32楼

xuelida

发表于 2012-4-12 19:13:23 |只看作者 |坛友微信交流群

这个程序对应的文章可以看我前面上传的一个pdf文件，还有就是Perron and Ng (1996), Review of Economic Studies, Useful
Modifications to Unit root tests with dependent errors and their local
asymptotic properties, vol 63, p.435-464.

proc(15)=mic1(y,p,penalty,kmax,kmin);    @调用主程序@
local nt,yt,ahat,r,fit,s2u,sumyt2,ta,sar,bt,za,mza,msb;
local adf,krule,sarf,ktild,ssra,ssr1,ytf,pt,mpt,mzt,mza2,adfols,a1ols;
local z,cbar,a1,s2ar1,mza1,za1,yols,yols2,ahatols,rols,fitols,s2uols,krule1;
nt=rows(y);
if p == 0;             @p: 回归多项式，0对应着只有常数项，1对应常数项和线性趋势项@
z=ones(nt,1);
cbar=-7.0;
endif;
if p == 1;
z=ones(nt,1)~seqa(1,1,nt);
cbar=-13.5;
endif;
cbar=cbar;
{yt,ssra}=glsd(y,z,cbar);                               /* transforming the data*/
{ahat,r,fit}=olsqr2(yt[2:nt,1],yt[1:nt-1,1]); /* estimate the alpha-hat*/
s2u=r'r/(nt-1);
sumyt2=sumc(yt[1:nt-1,1]^2)/(nt-1)^2;
krule=s2ar(yt,penalty,kmax,kmin);                /* estimate s2ar */
{adf,a1,sar}=adfp(yt,krule);
/*constructing zalpha and the M tests using GLS detrended data*/
bt=nt-1;
za=bt*(ahat-1)-(sar-s2u)/(2*sumyt2);
mza=((yt[nt,1]^2)/bt-sar)/(2*sumyt2);
msb=sqrt(sumyt2/sar);
mzt=mza*msb;
/* construct PT and MPT */
{ytf,ssr1}=glsd(y,z,0);
pt=(ssra-(1+cbar/nt)*ssr1)/sar;
if p == 0;
mpt=(cbar*cbar*sumyt2-cbar*(yt[nt,1]^2)/nt)/sar;
endif;
if  p == 1;
mpt=(cbar*cbar*sumyt2+(1-cbar)*(yt[nt,1]^2)/nt)/sar;
endif;

/* the ols detrending statistics */
{yols,yols2}=olsd(y,p);
{ahatols,rols,fitols}=olsqr2(yols[2:nt,1],yols[1:nt-1,1]);
s2uols=rols'rols/(nt-1);
krule1=s2ar(yols,penalty,kmax,kmin);
{adfols,a1ols,s2ar1}=adfp(yols,krule1);
za1=bt*(ahatols-1)-(s2ar1-s2uols)/(2*yols2);
mza1=((yols[nt,1]^2)/bt-s2ar1)/(2*yols2);
mza2=((yt[nt,1]^2)/bt-s2ar1)/(2*sumyt2);
if p == 0;
za = za| -8.35;
za1 = za1| -14.1;
mzt = mzt| -1.98;
mza = mza|-8.1;
mza1= mza1|-14.1;
mza2= mza2|-8.1;
msb = msb| .233;
adf = adf| -1.98;
mpt = mpt| 3.17;
pt = pt| 3.17;
adfols=adfols|-2.86;
else;
za = za| -17.3;
za1 = za1| -21.0;
mzt = mzt| -2.91;
mza1 = mza1|-21.3;
mza2 = mza2|-17.3;
mza = mza|-17.3;
msb = msb| .168;
adf = adf| -2.91;
mpt = mpt| 5.48;
pt = pt| 5.48;
adfols=adfols|-3.41;
endif;
retp(za,mza,msb,adf,pt,mpt,mzt,za1,mza1,mza2,adfols,krule,krule1,a1,ahatols);
endp;
proc(3)=adfp(yt,kstar);
local reg,dyt,i,rho,ee,ff,nef,s2e,xx,sre,adf,sumb,s2vec;
reg=lagn(yt,1);
dyt=diff(yt,1);
i=1;
do while i <= kstar;
reg=reg~lagn(dyt,i);
i=i+1;
endo;
dyt=trimr(dyt,kstar+1,0);
reg=trimr(reg,kstar+1,0);
{rho,ee,ff}=olsqr2(dyt,reg);
/*nef=rows(ee);*/
nef=rows(dyt);
s2e=ee'ee/nef;
xx=inv(reg'reg);
sre=xx[1,1]*s2e;
adf=rho[1,1]/sqrt(sre);
if kstar> 0;
sumb=sumc(rho[2:kstar+1]);
else;
sumb=0;
endif;
s2vec=s2e/(1-sumb)^2;
retp(adf,rho[1,1]+1,s2vec);
endp;

proc(1)=s2ar(yts,penalty,kmax,kmin);
local nt,min,s2vec,dyts,reg,k,b,e,fit,nef,s2e,dyts0,reg0,i;
local bic,kbic,sumb,j,msbar,gap,kopt,ssr,trgff,sumy,tau,mic,kk;
nt=rows(yts);
min=9999999999;
tau=zeros(kmax+1,1);
s2e=999*ones(kmax+1,1);

dyts=diff(yts,1);
reg=lagn(yts,1);
i=1;
do while i <= kmax;
reg=reg~lagn(dyts,i);
i=i+1;
endo;
dyts0=dyts;
reg0=reg;
/*loop over k*/
dyts0=trimr(dyts,kmax+1,0);
reg0=trimr(reg,kmax+1,0);
sumy=sumc(reg0[.,1].*reg0[.,1]);
nef=nt-kmax-1;
k=kmin;
do while k <= kmax;
/*{b,e,fit}=olsqr2(dyts0,reg0[.,1:k+1]);*/
b=dyts0/reg0[.,1:k+1];
e=dyts0-reg0[.,1:k+1]*b;
s2e[k+1]=e'e/nef;
tau[k+1]=(b[1]*b[1])*sumy/s2e[k+1];
k=k+1;
endo;
kk=seqa(0,1,kmax+1);
if penalty == 0;
mic=ln(s2e)+2.0*(kk+tau)./nef;
else;
mic=ln(s2e)+ln(nef)*(kk)./nef;
endif;
kopt=minindc(mic)-1;
retp(kopt);
endp;

proc(2)=glsd(y,z,cbar);
local nt,abar,ya,za,bhat,yt,ssr;
nt=rows(y);
abar=1+cbar/nt;
ya=zeros(nt,1);
za=zeros(nt,cols(z));
ya[1:1,1]=y[1:1,1];
za[1:1,.]=z[1:1,.];
ya[2:nt,1]=y[2:nt,1]-abar*y[1:nt-1,1];
za[2:nt,.]=z[2:nt,.]-abar*z[1:nt-1,.];
/*constructing the gls detrended series*/
bhat=inv(za'za)*za'ya;
yt=y-z*bhat;
ssr=(ya-za*bhat)'(ya-za*bhat);
retp(yt,ssr);
endp;
proc(1)= diff(x,k) ;
   if ( k == 0) ;
      retp(x) ;
   endif ;
   retp(zeros(k,cols(x))|trimr(x,k,0)-trimr(lagn(x,k),k,0)) ;
endp ;
proc(1)= lagn(x,n);
local y;
if n > 0;
y = zeros(n,cols(x))|trimr(x,0,n);
else;
y = trimr(x,abs(n),0)|zeros(abs(n),cols(x));
endif;
retp(y);
endp;
proc(2)=olsd(y,p);
local yd,reg,trend,ahat,fit,i;
reg=ones(rows(y),1);
trend=seqa(1,1,rows(y));
i=1;
do while i<=p;
reg=reg~trend^i;
i=i+1;
endo;
{ahat,yd,fit}=olsqr2(y,reg);
retp(yd,sumc(yd[1:rows(yd)-1]^2)/(rows(yd)-1)^2);
endp;

使用道具举报