Gauss 编程应用于学科前沿一：单位根

下面我们看看Phillips-Perron检验，考虑序列相关和异方差下的Newey-West估计量。
开始设置 e'e/n1 ，然后修正为2*wi*(e_i'e)/n1。
nw=e'e/n1;
i=1;
do while i<=lags;
    w=1-i/(lags+1);
    nw=nw+2*w*e[1:n1-i]'e[i+1:n1]/n1;
    i=i+1;
endo;
还是调用前面的数据

data=reshape(data,rows(data)/2,2);
y=data[.,1];
lags=5;
n=rows(y);
dy=y[2:n]-y[1:n-1];
x=y[1:n-1];
n1=rows(x);
b=inv(x'x)*x'dy;
e=dy-x*b;
nw=e'e/n1;
i=1;
do while i<=lags;
   w=1-i/(lags+1);
   nw=nw+2*w*e[1:n1-i]'e[i+1:n1]/n1;
   i=i+1;
endo;
print "NW estimator=" nw;

用这个程序作为长期方差的nw估计，则DF统计量修正为，
taupp=sqrt(s2u/s2tq)*tau-1/2*(s2tq-s2u)/sqrt(s2tq)*n1*setau/sig其中，
s2u=e'e/n1 ， s2tq=nw
setau 是第一个系数的标准误，sig是残差向量e的标准偏。

data=reshape(data,rows(data)/2,2);
y=data[.,1];
lags=5;
call pp(y,lags);

proc(0)=pp(y,lags);
   local n,dy,x,n1,k,b,e,s2hat,varb,se,tau,setau,sig,nw,w,i,taupp,s2u,s2tq;
   n=rows(y);
   dy=y[2:n]-y[1:n-1];
   x=y[1:n-1];
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1]; print "tau  =" tau;
   setau=se[1];
   sig=stdc(e);
   nw=e'e/n1;
   i=1;
   do while i<=lags;
      w=1-i/(lags+1);
      nw=nw+2*w*e[1:n1-i]'e[i+1:n1]/n1;
      i=i+1;
   endo;
   s2u=e'e/n1;
   s2tq=nw;
   taupp=sqrt(s2u/s2tq)*tau-1/2*(s2tq-s2u)/sqrt(s2tq)*n1*setau/sig;
   print "taupp=" taupp;;print/rz "( n=" n1 ")";
endp;

如果滞后长度=0，则tau和taupp相同。

下面把我们前面选择滞后长度的拇指准则（rule of thumb ）用在pp检验
模型1：不带趋势项和漂移项    x=y[1:n-1];
模型2：不带趋势项                  x=y[1:n-1]~ones(n-1,1);
模型3：都带                            x=y[1:n-1]~ones(n-1,1)~seqa(1,1,n-1);

dy=y[2:n]-y[1:n-1];


data=reshape(data,rows(data)/2,2);
y=data[.,1];
print "Model I:";
call pp1(y);
print "Model II(with drift):";
call pp2(y);
print "Model III(with drift and trend):";
call pp3(y);

proc(0)=pp1(y);
   local n,lags,dy,x,n1,k,b,e,s2hat,varb,se,tau,setau,sig,nw,w,i,taupp,s2u,s2tq;
   n=rows(y);
   lags=floor(4*(n/100)^(1/4));
   dy=y[2:n]-y[1:n-1];
   x=y[1:n-1];
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1]; print "tau  =" tau;
   setau=se[1];
   sig=stdc(e);
   nw=e'e/n1;
   i=1;
   do while i<=lags;
      w=1-i/(lags+1);
      nw=nw+2*w*e[1:n1-i]'e[i+1:n1]/n1;
      i=i+1;
   endo;
   s2u=e'e/n1;
   s2tq=nw;
   taupp=sqrt(s2u/s2tq)*tau-1/2*(s2tq-s2u)/sqrt(s2tq)*n1*setau/sig;
   print "taupp=" taupp;;print/rz "( n=" n1 ")";
endp;

proc(0)=pp2(y);
   local n,lags,dy,x,n1,k,b,e,s2hat,varb,se,tau,setau,sig,nw,w,i,taupp,s2u,s2tq;
   n=rows(y);
   lags=floor(4*(n/100)^(1/4));
   dy=y[2:n]-y[1:n-1];
   x=y[1:n-1]~ones(n-1,1);
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1]; print "tau  =" tau;
   setau=se[1];
   sig=stdc(e);
   nw=e'e/n1;
   i=1;
   do while i<=lags;
      w=1-i/(lags+1);
      nw=nw+2*w*e[1:n1-i]'e[i+1:n1]/n1;
      i=i+1;
   endo;
   s2u=e'e/n1;
   s2tq=nw;
   taupp=sqrt(s2u/s2tq)*tau-1/2*(s2tq-s2u)/sqrt(s2tq)*n1*setau/sig;
   print "taupp=" taupp;;print/rz "( n=" n1 ")";
endp;

proc(0)=pp3(y);
   local n,lags,dy,x,n1,k,b,e,s2hat,varb,se,tau,setau,sig,nw,w,i,taupp,s2u,s2tq;
   n=rows(y);
   lags=floor(4*(n/100)^(1/4));
   dy=y[2:n]-y[1:n-1];
   x=y[1:n-1]~ones(n-1,1)~seqa(1,1,n-1);
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1]; print "tau  =" tau;
   setau=se[1];
   sig=stdc(e);
   nw=e'e/n1;
   i=1;
   do while i<=lags;
      w=1-i/(lags+1);
      nw=nw+2*w*e[1:n1-i]'e[i+1:n1]/n1;
      i=i+1;
   endo;
   s2u=e'e/n1;
   s2tq=nw;
   taupp=sqrt(s2u/s2tq)*tau-1/2*(s2tq-s2u)/sqrt(s2tq)*n1*setau/sig;
   print "taupp=" taupp;;print/rz "( n=" n1 ")";
endp;

下面开始进行DF-GLS去势检验

/*DF-GLS 估计，利用rho=1-7/n去除原序列的趋势 */
new; cls;
let data=
              80         10236000
              78         10351000
              71          9371000
              73          9126000
              70          8230000
              71         12111000
              72          5364000
              76          6222000
              76          3952000
              78          3859000
              80          4798000
              80          5133000
              77          4240000
              69         13687000
              69          5248000
              68          4058000
              67          5253000
              66          8963000
              67          6457000
              70          8882000
              69          4134000
              68          3736000
              69          3887000
              69          2095000
              68          1980000
              63          7433000
              63          5487000
              63          3996000
              58         10858000
              60         12095000
              57          7251000
              53         17632000
              55         42481000
              65         16051000
              64          9145000
              62         11264000
              64          4799000
              64          6888000
              66          7093000
              69          6598000
              70          6280000
              69          5372000
              65          7381000
              64          2792000
              63          4258000
              60          3074000
              62          4427000
              60          3767000
              61          3041000
              60         10240000
;
data=reshape(data,rows(data)/2,2);
y=data[.,1];
n=rows(y);
dy=y[2:n]-y[1:n-1];
y_1=y[1:n-1];
print "DF test statistic:";
call dftau(dy,y_1~ones(n-1,1));
print "DF-GLS test statistic:";
call dfgls(y);

proc dfgls(y);
local n,rho,y1,x,b,dy,y_1,tau;
   n=rows(y);
   rho=1-7/n;                        @去势因子@
   y1=y[1]|(y[2:n]-rho*y[1:n-1]);    @构造去势的y@
   x=1|(1-rho)*ones(n-1,1);          @构造x{1,1-rho,1-rho,...}@
   b=inv(x'x)*x'y1;                  @利用y和x回归得到回归系数b@
   y=y-x*b;                          
   dy=y[2:n]-y[1:n-1];               @差分@
   y_1=y[1:n-1];                     @y的一阶滞后@
   tau=dftau(dy,y_1);
   retp(tau);
endp;
proc dftau(dy,x);
local n1,k,b,e,s2hat,varb,se,tau;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1];
   print "tau=" tau;;
   print/rz "( n=" n1 ")";
   retp(tau);
endp;

楼主非常赞，支持！

上面我们虽然得到了检验统计量的值，但是我们没有得到临界值，DF-GLS临界值可以由下面得到
new; cls;
pp={0.10,0.05,0.025,0.01};
times=10000;                     /* increase more in normal version of GAUSS */
label={"N" "10%" "5%" "2.5%" "1%"};
print $label;
print/rz 49~dfglssta(49,pp,times)';

proc dfglssta(n,pp,times);
local tau,j,e,y,stat;
   tau=zeros(times,1);
   j=1;
   do while j<=times;
      e=rndn(n,1);
      y=cumsumc(e);
      tau[j]=dfgls(y);
      j=j+1;
   endo;
   tau=sortc(tau,1);                  @排序@
   stat=tau[round(pp*times)];   @百分位数@
   retp(stat);
endp;
proc dfgls(y);
   local n,rho,y1,x,b,dy,y_1,tau;
   n=rows(y);
   rho=1-7/n;
   y1=y[1]|(y[2:n]-rho*y[1:n-1]);
   x=1|(1-rho)*ones(n-1,1);
   b=inv(x'x)*x'y1;
   y=y-x*b;
   dy=y[2:n]-y[1:n-1];
   y_1=y[1:n-1];
   tau=dftau(dy,y_1);
   retp(tau);
endp;
proc dftau(dy,x);
   local n1,k,b,e,s2hat,varb,se,tau;
   n1=rows(x);
   k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1];
   retp(tau);
endp;

上面程序模拟的百分位数结果为

              N          10%            5%         2.5%           1%
              49       -1.652        -1.976       -2.274       -2.629

对应的就可以当做临界值使用

如果带有时间趋势项，根据perron的文章可以使用rho=1-13.5/n，可以构造解释变量矩阵
   x1=1|(1-rho)*ones(n-1,1);
   xt=t[1]|(t[2:n]-rho*t[1:n-1]);
   x=x1~xt;

DF检验带有drift and trend 项，数据见前例
读者可以看看perron的文章，在附件里

new; cls;
data=reshape(data,rows(data)/2,2);
y=data[.,1];

n=rows(y);
dy=y[2:n]-y[1:n-1];
y_1=y[1:n-1];
print "DF test statistic:";
call dftau(dy,y_1~ones(n-1,1)~seqa(1,1,n-1));

print "DF-GLS test statistic:";
call dfglst(y);


proc dfglst(y);
   local n,rho,t,y1,x1,xt,x,b,dy,y_1,tau;
   n=rows(y);
   rho=1-13.5/n;
   t=seqa(1,1,n);
   y1=y[1]|(y[2:n]-rho*y[1:n-1]);
   x1=1|(1-rho)*ones(n-1,1);
   xt=t[1]|(t[2:n]-rho*t[1:n-1]);
   x=x1~xt;
   b=inv(x'x)*x'y1;
   y=y-x*b;
   dy=y[2:n]-y[1:n-1];
   y_1=y[1:n-1];
   tau=dftau(dy,y_1);
   retp(tau);
endp;

proc dftau(dy,x);
   local n1,k,b,e,s2hat,varb,se,tau;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   tau=b[1]./se[1];
   print "tau=" tau;; print/rz "( n=" n1 ")";
   retp(tau);
endp;

读者可以利用模拟百分位数的程序，改成带drift and trend 项，重新计算临界值，做检验

ADF-GLS检验，使用MAIC 准则选择滞后长度
数据见前例，下面的例子没有趋势项t，读者可以照着前面的例子，加入趋势项t

new; cls;
data=reshape(data,rows(data)/2,2);
y=data[.,1];

print "ADF test:";
call adflags(y,10);   @设置最大滞后值10@
print;
print "ADF-GLS(level) test:";
call adfglslags(y,10);


proc(0)=adflags(y,maxlags);        @ ADF test  @      
local n,dy0,y_1,dy,n1,k,x,tau,dy_i,i,label,b,e,sig2p,taup,maic;
   maic=miss(zeros(maxlags+1,1),0);
   n=rows(y);
   dy0=y[2:n]-y[1:n-1];
   y_1=y[1:n-1];
   /* lag 0 */
   dy=dy0[maxlags+1:n-1];
   y_1=y_1[maxlags+1:n-1];
   n1=rows(y_1);
   k=2;
   x=y_1~ones(n1,1);
   tau=tOLS(dy,x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   sig2p=e'e/n1;
   taup=(b[1]^2)*(y_1'y_1)/sig2p;
   maic[1]=ln(sig2p)+2*(taup+0)/n1;         
   print/lz "n=" n1;
   label={"lags" "tau" "t dy_max" "p-value" "MAIC"};
   print $label;
   print/rz 0~tau[1]~miss(0,0)~miss(0,0)~maic[1];
   /* lag 1 to lags */
   dy_i=ones(n1,1);
   i=1;
   do while i<=maxlags;
      dy_i=dy_i~dy0[maxlags+1-i:n-1-i];
      x=y_1~dy_i;
      k=cols(x);
      tau=tOLS(dy,x);
      b=inv(x'x)*x'dy;
      e=dy-x*b;
      sig2p=e'e/n1;
      taup=(b[1]^2)*(y_1'y_1)/sig2p;
      maic[i+1]=ln(sig2p)+2*(taup+i)/n1;
      print/rz i~tau[1]~tau[k]~2*cdftc(abs(tau[k]),n1-k)~maic[i+1];
      i=i+1;
   endo;
   print/lz "# of opt.lags=" minindc(maic)-1;
endp;

proc(0)=adfglslags(y,maxlags);         @ ADF-GLS(level) test@
   local n,rho,y1,dy0,y_1,dy,n1,k,x,tau,dy_i,i,label,b,e,sig2p,taup,maic;
   maic=miss(zeros(maxlags+1,1),0);
   n=rows(y);
   /* detrending */
   rho=1-7/n;
   y1=y[1]|(y[2:n]-rho*y[1:n-1]);
   x=1|(1-rho)*ones(n-1,1);
   b=inv(x'x)*x'y1;
   y=y-x*b;
   dy0=y[2:n]-y[1:n-1];
   y_1=y[1:n-1];
   /* lag 0 */
   dy=dy0[maxlags+1:n-1];
   y_1=y_1[maxlags+1:n-1];
   n1=rows(y_1);
   k=1;
   x=y_1;
   tau=tOLS(dy,x);
   b=inv(x'x)*x'dy;
   e=dy-x*b;
   sig2p=e'e/n1;
   taup=(b[1]^2)*(y_1'y_1)/sig2p;
   maic[1]=ln(sig2p)+2*(taup+0)/n1;
   print/lz "n=" n1;
   label={"lags" "tau" "t dy_max" "p-value" "MAIC"};
   print $label;
   print/rz 0~tau[1]~miss(0,0)~miss(0,0)~maic[1];
   /* lag 1 to lags */
   i=1;
      dy_i=dy0[maxlags+1-i:n-1-i];
      x=y_1~dy_i;
      k=cols(x);
      tau=tOLS(dy,x);
      b=inv(x'x)*x'dy;
      e=dy-x*b;
      sig2p=e'e/n1;
      taup=(b[1]^2)*(y_1'y_1)/sig2p;
      maic[i+1]=ln(sig2p)+2*(taup+i)/n1;
      print/rz i~tau[1]~tau[k]~2*cdftc(abs(tau[k]),n1-k)~maic[i+1];
   i=2;
   do while i<=maxlags;
      dy_i=dy_i~dy0[maxlags+1-i:n-1-i];
      x=y_1~dy_i;
      k=cols(x);
      tau=tOLS(dy,x);
      b=inv(x'x)*x'dy;
      e=dy-x*b;
      sig2p=e'e/n1;
      taup=(b[1]^2)*(y_1'y_1)/sig2p;
      maic[i+1]=ln(sig2p)+2*(taup+i)/n1;
      print/rz i~tau[1]~tau[k]~2*cdftc(abs(tau[k]),n1-k)~maic[i+1];
      i=i+1;
   endo;
   print/lz "# of opt.lags=" minindc(maic)-1;
endp;

proc tOLS(y,x);
   local n1,k,b,e,s2hat,varb,se,t;
   n1=rows(x); k=cols(x);
   b=inv(x'x)*x'y;
   e=y-x*b;
   s2hat=e'e/(n1-k);
   varb=s2hat*inv(x'x);
   se=sqrt(diag(varb));
   t=b./se;
   retp(t);
endp;

上面程序的结果为
ADF test:
n= 39              
            lags              tau         t dy_max          p-value             MAIC
               0       -2.6544716                .                .        52.653032
               1       -2.7598758       0.91812825       0.36466239        59.725523
               2       -2.4888546       -0.6546252       0.51698824        53.509516
               3        -2.602916       0.91707459       0.36556321        61.745706
               4       -2.5715032       0.42385468       0.67442197        69.001319
               5       -2.4329103     0.0065598943       0.99480673        69.153789
               6       -2.2166857      -0.26600835       0.79199404         64.90662
               7       -2.2379659       0.47418382       0.63880209          70.3685
               8       -1.9263461       -1.2376599       0.22576958        57.165185
               9       -1.9301568      -0.98428148       0.33340647        59.371316
              10       -1.9362932      -0.63686844        0.5295738        62.021903
# of opt.lags= 0               

ADF-GLS(level) test:
n= 39              
            lags              tau         t dy_max          p-value             MAIC
               0       -1.3545972                .                .        2.1969967
               1        -1.270363       0.44158716       0.66135806        2.2352195
               2       -1.3949018       -1.2349445       0.22485086        2.2673229
               3       -1.3169155       0.44227396       0.66101204        2.3052449
               4        -1.321108      -0.37378498        0.7108851        2.3565465
               5       -1.3350658      -0.63935143       0.52700796        2.4013202
               6       -1.3632774      -0.88060291        0.3850989        2.4373618
               7       -1.3268783      0.080486454       0.93636775        2.4869745
               8       -1.4055136        -1.628072       0.11396808        2.4729003
               9       -1.5507825      -0.96533816       0.34235713        2.5343324
              10       -1.5639181       0.37785183          0.70839        2.5935194
# of opt.lags= 0