{"word_count": 1612, "character_count": 13500}
This week, in the MAT8181 Time Series course, we've discussed unit root tests. According to Wold's theorem, if  is  (weakly) stationnary thenwhere  is the innovation process, and where  is some deterministic series (just to get a result as general as possible). Observe thatas discussed in a previous post. To go one step further, there is also the Beveridge-Nelson decomposition : an integrated of order one process, defined ascan be represented asa linear trend  a random walk  a stationary remaining termi.e.where  is the polynomial with terms , whereFor unit-root tests, we will use various representation of the process. In order to illustrate the implementation of those tests, consider the following series&#62; E=rnorm(240)&#62; X=cumsum(E)&#62; plot(X,type="l")Dickey Fuller (standard)Here, for the simple version of the Dickey-Fuller test, we assume thatand we would like to test if  (or not). We can write the previous representation asso we simply have to test if the regression coefficient in the linear regression is - or not - null. Which can be done with Student's test. If we consider the previous model without the linear drift, we have to consider the following regression&#62; lags=0&#62; z=diff(X)&#62; n=length(z)&#62; z.diff=embed(z, lags+1)[,1]&#62; z.lag.1=X[(lags+1):n]&#62; summary(lm(z.diff~0+z.lag.1 ))Call:lm(formula = z.diff ~ 0 + z.lag.1)Residuals:     Min       1Q   Median       3Q      Max -2.84466 -0.55723 -0.00494  0.63816  2.54352 Coefficients:         Estimate Std. Error t value Pr(&#62;|t|)z.lag.1 -0.005609   0.007319  -0.766    0.444Residual standard error: 0.963 on 238 degrees of freedomMultiple R-squared:  0.002461,Adjusted R-squared:  -0.00173 F-statistic: 0.5873 on 1 and 238 DF,  p-value: 0.4442Our testing procedure will be based on the Student's t value,&#62; summary(lm(z.diff~0+z.lag.1 ))$coefficients[1,3][1] -0.7663308which is exactly the value computed using&#62; library(urca)&#62; df=ur.df(X,type="none",lags=0)&#62; df############################################################### # Augmented Dickey-Fuller Test Unit Root / Cointegration Test # ############################################################### The value of the test statistic is: -0.7663The interpretation of this value can be done using critical values (99%, 95%, 90%)&#62; qnorm(c(.01,.05,.1)/2)[1] -2.575829 -1.959964 -1.644854If the statistics exceeds those values, then the series is not stationnary, since we cannot reject the assumption that . So we might conclude that there is a unit root. Actually, those critical values are obtained using&#62; summary(df)############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression none Call:lm(formula = z.diff ~ z.lag.1 - 1)Residuals:     Min       1Q   Median       3Q      Max -2.84466 -0.55723 -0.00494  0.63816  2.54352 Coefficients:         Estimate Std. Error t value Pr(&#62;|t|)z.lag.1 -0.005609   0.007319  -0.766    0.444Residual standard error: 0.963 on 238 degrees of freedomMultiple R-squared:  0.002461,Adjusted R-squared:  -0.00173 F-statistic: 0.5873 on 1 and 238 DF,  p-value: 0.4442Value of test-statistic is: -0.7663 Critical values for test statistics:       1pct  5pct 10pcttau1 -2.58 -1.95 -1.62The problem with R is that there are several packages that can be used for unit root tests. Just to mention another one,&#62; library(tseries)&#62; adf.test(X,k=0)Augmented Dickey-Fuller Testdata:  XDickey-Fuller = -2.0433, Lag order = 0, p-value = 0.5576alternative hypothesis: stationaryWe do have here also a test where the null hypothesis is that there is a unit root. But the p-value is quite different. What is odd is that we have&#62; 1-adf.test(X,k=0)$p.value[1] 0.4423705&#62; df@testreg$coefficients[4][1] 0.4442389(but I think it is a coincidence).Augmented Dickey FullerIt is possible to had some lags in the regression. For instance, we can considerAgain, we have to check if one coefficient is null, or not. And this can be done using Student's t test.&#62; lags=1&#62; z=diff(X)&#62; n=length(z)&#62; z.diff=embed(z, lags+1)[,1]&#62; z.lag.1=X[(lags+1):n]&#62; k=lags+1&#62; z.diff.lag = embed(z, lags+1)[, 2:k]&#62; summary(lm(z.diff~0+z.lag.1+z.diff.lag ))Call:lm(formula = z.diff ~ 0 + z.lag.1 + z.diff.lag)Residuals:     Min       1Q   Median       3Q      Max -2.87492 -0.53977 -0.00688  0.64481  2.47556 Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)z.lag.1    -0.005394   0.007361  -0.733    0.464z.diff.lag -0.028972   0.065113  -0.445    0.657Residual standard error: 0.9666 on 236 degrees of freedomMultiple R-squared:  0.003292,Adjusted R-squared:  -0.005155 F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777&#62; summary(lm(z.diff~0+z.lag.1+z.diff.lag ))$coefficients[1,3][1] -0.7328138This value is the one obtained using&#62; df=ur.df(X,type="none",lags=1)&#62; summary(df)############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression none Call:lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)Residuals:     Min       1Q   Median       3Q      Max -2.87492 -0.53977 -0.00688  0.64481  2.47556 Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)z.lag.1    -0.005394   0.007361  -0.733    0.464z.diff.lag -0.028972   0.065113  -0.445    0.657Residual standard error: 0.9666 on 236 degrees of freedomMultiple R-squared:  0.003292,Adjusted R-squared:  -0.005155 F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777Value of test-statistic is: -0.7328 Critical values for test statistics:       1pct  5pct 10pcttau1 -2.58 -1.95 -1.62And again, other pckages can be used:&#62; adf.test(X,k=1)Augmented Dickey-Fuller Testdata:  XDickey-Fuller = -1.9828, Lag order = 1, p-value = 0.5831alternative hypothesis: stationaryHopefully, the conclusion is the same (we should reject the assumption that the series is stationary, but I am not sure about the computation of the p-value).Augmented Dickey Fuller with trend and driftSo far, we have not included the drift in our model. But this is simple to do (this will be called the augmented version of the previous procedure): we just have to include a constant in the regression,&#62; summary(lm(z.diff~1+z.lag.1+z.diff.lag ))Call:lm(formula = z.diff ~ 1 + z.lag.1 + z.diff.lag)Residuals:     Min       1Q   Median       3Q      Max -2.91930 -0.56731 -0.00548  0.62932  2.45178 Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)  (Intercept)  0.29175    0.13153   2.218   0.0275 *z.lag.1     -0.03559    0.01545  -2.304   0.0221 *z.diff.lag  -0.01976    0.06471  -0.305   0.7603  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.9586 on 235 degrees of freedomMultiple R-squared:  0.02313,Adjusted R-squared:  0.01482 F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393The statistics of interest are obtained here considering some analysis of variance outputs, where this model is compared with the one without the integrated part, and the drift,&#62; summary(lm(z.diff~1+z.lag.1+z.diff.lag ))$coefficients[2,3][1] -2.303948&#62; anova(lm(z.diff ~ z.lag.1 + 1 + z.diff.lag),lm(z.diff ~ 0 + z.diff.lag))$F[2][1] 2.732912Those two values are the ones obtained also with&#62; df=ur.df(X,type="drift",lags=1)&#62; summary(df)############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression drift Call:lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)Residuals:     Min       1Q   Median       3Q      Max -2.91930 -0.56731 -0.00548  0.62932  2.45178 Coefficients:            Estimate Std. Error t value Pr(&#62;|t|)  (Intercept)  0.29175    0.13153   2.218   0.0275 *z.lag.1     -0.03559    0.01545  -2.304   0.0221 *z.diff.lag  -0.01976    0.06471  -0.305   0.7603  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.9586 on 235 degrees of freedomMultiple R-squared:  0.02313,Adjusted R-squared:  0.01482 F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393Value of test-statistic is: -2.3039 2.7329 Critical values for test statistics:       1pct  5pct 10pcttau2 -3.46 -2.88 -2.57phi1  6.52  4.63  3.81And we can also include a linear trend,&#62; temps=(lags+1):n&#62; summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))Call:lm(formula = z.diff ~ 1 + temps + z.lag.1 + z.diff.lag)Residuals:     Min       1Q   Median       3Q      Max -2.87727 -0.58802 -0.00175  0.60359  2.47789 Coefficients:              Estimate Std. Error t value Pr(&#62;|t|)  (Intercept)  0.3227245  0.1502083   2.149   0.0327 *temps       -0.0004194  0.0009767  -0.429   0.6680  z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.9603 on 234 degrees of freedomMultiple R-squared:  0.0239,Adjusted R-squared:  0.01139 F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287&#62; summary(lm(z.diff~1+temps+z.lag.1+z.diff.lag ))$coefficients[3,3][1] -1.98282&#62; anova(lm(z.diff ~ z.lag.1 + 1 + temps+ z.diff.lag),lm(z.diff ~ 1+ z.diff.lag))$F[2][1] 2.737086while R function returns&#62; df=ur.df(X,type="trend",lags=1)&#62; summary(df)############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression trend Call:lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)Residuals:     Min       1Q   Median       3Q      Max -2.87727 -0.58802 -0.00175  0.60359  2.47789 Coefficients:              Estimate Std. Error t value Pr(&#62;|t|)  (Intercept)  0.3227245  0.1502083   2.149   0.0327 *z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *tt          -0.0004194  0.0009767  -0.429   0.6680  z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243  ---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.9603 on 234 degrees of freedomMultiple R-squared:  0.0239,Adjusted R-squared:  0.01139 F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287Value of test-statistic is: -1.9828 1.8771 2.7371 Critical values for test statistics:       1pct  5pct 10pcttau3 -3.99 -3.43 -3.13phi2  6.22  4.75  4.07phi3  8.43  6.49  5.47KPSS testHere, in the KPSS testing procedure, two models can be considerd : with a drift, or with a linear trend. Here, the null hypothesis is that the series is stationnary.With a drift, the code is&#62; summary(ur.kpss(X,type="mu"))####################### # KPSS Unit Root Test # ####################### Test is of type: mu with 4 lags. Value of test-statistic is: 0.972 Critical value for a significance level of:                 10pct  5pct 2.5pct  1pctcritical values 0.347 0.463  0.574 0.73while it will be, in the case there is a trend&#62; summary(ur.kpss(X,type="tau"))####################### # KPSS Unit Root Test # ####################### Test is of type: tau with 4 lags. Value of test-statistic is: 0.5057 Critical value for a significance level of:                 10pct  5pct 2.5pct  1pctcritical values 0.119 0.146  0.176 0.216One more time, it is possible to use another package to get the same test (but again, a different output)&#62; kpss.test(X,"Level")KPSS Test for Level Stationaritydata:  XKPSS Level = 1.1997, Truncation lag parameter = 3, p-value = 0.01&#62; kpss.test(X,"Trend")KPSS Test for Trend Stationaritydata:  XKPSS Trend = 0.6234, Truncation lag parameter = 3, p-value = 0.01At least, there is some kind of consistency, since we keep rejecting the stationnary assumption, for that series.Philipps-Perron testThe Philipps-Perron test is based on the ADF procedure. The code is here&#62; PP.test(X)Phillips-Perron Unit Root Testdata:  XDickey-Fuller = -2.0116, Truncation lag parameter = 4, p-value = 0.571with again, a possible alternative with the other package&#62; pp.test(X)Phillips-Perron Unit Root Testdata:  XDickey-Fuller Z(alpha) = -7.7345, Truncation lag parameter = 4, p-value= 0.6757alternative hypothesis: stationary ComparisonI will not spend more time comparing the different codes, in R, to run those tests. Let us spend some additional time on a quick comparison of those three procedure. Let us generate some autoregressive processes, with more or less autocorrelation, as well as some random walk, and let us see how those tests perform :&#62; n=100&#62; AR=seq(1,.7,by=-.01)&#62; P=matrix(NA,3,31)&#62; M1=matrix(NA,1000,length(AR))&#62; M2=matrix(NA,1000,length(AR))&#62; M3=matrix(NA,1000,length(AR))&#62; for(i in 1:(length(AR)+1)){+ for(s in 1:1000){+ if(i==1) X=cumsum(rnorm(n))+ if(i!=1) X=arima.sim(n=n,list(ar=AR[i]))+ library(urca)+ M2[s,i]=as.numeric(pp.test(X)$p.value)+ M1[s,i]=as.numeric(kpss.test(X)$p.value)+ M3[s,i]=as.numeric(adf.test(X)$p.value)+ }}Here, we would like to count how many times the p-value of our tests exceed 5%,&#62; prop05=function(x) mean(x&#62;.05)+ P[1,]=1-apply(M1,2,prop05)+ P[2,]=apply(M2,2,prop05)+ P[3,]=apply(M3,2,prop05)+ }&#62; plot(AR,P[1,],type="l",col="red",ylim=c(0,1),ylab="proportion of non-stationnary + series",xlab="autocorrelation coefficient")&#62; lines(AR,P[2,],type="l",col="blue")&#62; lines(AR,P[3,],type="l",col="green")&#62; legend(.7,1,c("ADF","KPSS","PP"),col=c("green","red","blue"),lty=1,lwd=1)&#160;We can see here how poorly Dickey-Fuller test behave, since a 50% (at least) of our autoregressive processes are considered as non-stationnary.