function result = phillips(y,x,flag)% PURPOSE: compute Phillips-Perron test of the unit-root hypothesis%          based on a Dickey-Fuller/Augmented Dickey-Fuller regression% -------------------------------------------------------------------% USAGE: results = phillips(y,x,flag)%  where:              y = dependent variable (time series vector in levels)%                      x = explanatory variables matrix%                   flag = 0 for DF regression (default)%                        = nlag the # of lags for ADF regression % -------------------------------------------------------------------% RETURNS: a structure variable result%          results.meth = 'phillips'%          results.pstat = Phillips-Peron statistic, autocorrelation/heteroskedasticity %                         corrected t-ratio for the unit-root coefficient based on%                         a Dickey-Fuller or Augmented Dickey-Fuller regression%          results.crit  = (6 x 1) vector of critical values%                        [1% 5% 10% 90% 95% 99%] quintiles    %          results.nlag  = nlag, (should equal flag input)  %          results.alpha = autoregressive parameter estimate% -------------------------------------------------------------------% NOTES:    reject the null of a unit root, even in the presence of serial correlation%           and/or heteroskedasticity, if result.pstat is statistically significant.% -------------------------------------------------------------------% REFERENCES:  Phillips, Peter (1987), "Time Series Regression with a Unit Root", % Econometrica, vol. 55, no. 2 (March), pp. 277-301% Phillips, Peter & Pierre Perron (1988), "Testing for a Unit Root in Time Series%    Regression", Biometrika, vol. 75, no. 2 (June), pp. 335-346% -------------------------------------------------------------------% written by:% James P. LeSage, Dept of Economics% University of Toledo% 2801 W. Bancroft St,% Toledo, OH 43606% jlesage@spatia-econometrics.com% input checkingif nargin == 2, % use DF as defaultflag = 0;end;if nargin < 2 | nargin > 3error('Wrong # of arguments to phillips');end;result.meth = 'phillips';nobs = length(y); if flag == 0, % use DF     b = inv(x'*x)*x'*y;     r       = y - x*b;     dep     = tdiff(r,1);     dep = trimr(dep,1,0);     k       = 0     ;     z       = trimr(lag(r,1),1,0) ;   beta      = z\dep;   %resid       = dep - z*beta;     % BUG fix suggested by % Nick Firoozye% Sanford C. Bernstein, Inc% 767 Fifth Avenue, #21-49% New York, NY 10153     resid = detrend(dep,0)- detrend(z,0)*beta;      sigma      = (resid'*resid)/(rows(dep)-cols(z));     var_cov = sigma*inv(z'*z) ;else, % use ADF     b = inv(x'*x)*x'*y;     r       = y - x*b;     dep     = tdiff(r,1);     dep = trimr(dep,1,0);     k       = 1     ;     z       = trimr(lag(r,1),1,0) ;l = flag;          while (k <= l)           z = [z lag(dep,k)];           k = k + 1 ;          end;     z       = trimr(z,l,0) ;     dep     = trimr(dep,l,0) ;     beta    = detrend(z,0)\detrend(dep,0) ;     % resid     = dep - z*beta ;     % BUG fix suggested by % Nick Firoozye% Sanford C. Bernstein, Inc% 767 Fifth Avenue, #21-49% New York, NY 10153     resid = detrend(dep,0)- detrend(z,0)*beta;      sigma      = (resid'*resid)/(rows(dep)-cols(z));     var_cov = sigma*inv(z'*z) ;end;rss = resid'*resid;stdb = sqrt(var_cov(1,1));tstat = beta(1,1)/stdb;% The number of autocovariances of the residuals to be used, % based on ad-hoc rule of thumbcovs = fix(4*(nobs/100)^0.25);% Apply Newey & West (1987) adjustment:nw = rss;for j = 1 : covs   nw = nw + 2*(1-j/(covs+1))*(resid(1:end-j)'*resid(j+1:end));endnw = nw/nobs;% Obtain the Phillips-Perron statistic, which follows the% tabulated Dickey-Fuller distribution, and obtain its statistical significance:result.pstat = sqrt(rss/nobs / nw) * tstat - 1/2* (nw - rss/nobs)/sqrt(nw) * nobs*stdb/sigma;result.crit = ztcrit(nobs,0);result.alpha = beta(1,1);result.nlag = flag;