5.

STATIONARY STOCHASTIC PROCESSES

5.1 Definition
According to Gujarati and Porter (2009: p.740), a key concept underlying stochastic
process that has received a great attention and investigation by time series analysts is
the stationary stochastic process. Broadly speaking, “a stochastic process is said to be
stationary if its mean and variance are constant over time and the value of the
covariance24 between the two periods depends only on the distance or gap or lag
between the two time periods and not the actual time at which the covariance is
computed” (Gujarati and Porter, 2009: p.740)25. In the time series literature, such a
stochastic process is known as a weakly stationary or covariance stationary. By
contrast, a time series is strictly stationary26 if all the moments of its probability
distribution are time-invariant over time. If, however, the stationary process is normal,
the weakly stationary stochastic process is also strictly stationary (Gujarati and Porter,
2009: p.740). For most practical applications, the weak type of stationarity generally
suffices. According to Kocenda and Cerny (2017: p.17), the most frequently used
stationarity concept in econometrics is the concept of covariance stationarity. Therefore,
throughout this series of notes, we will for simplicity usually use only the term
stationarity instead of covariance or weak stationarity. According to Asteriou and Hall
(2016: p.277), a weakly stationary series is characterized by:
(a) exhibits mean reversion in that it fluctuates around a constant long-run mean;
(b) has a finite variance that is time-invariant; and
(c) has a theoretical correlogram27 that diminishes as the lag length increases.
In its simplest terms a time series Yt is said to be (weakly) stationary (hereafter refer to
stationary) if it has the following properties (Asteriou and Hall, 2016: p.277; Gujarati
and Porter, 2009: p.740; Hill et al., 2018: p.566; Kocenda and Cerny, 2017: p.17):
(a) mean: E(Yt) =  (constant for all t) (5.1)
i.e., t = t+k =  <  for all t, k
(b) variance: Var(Yt) = E(Yt – )2 = 2 (constant for all t) (5.2)
i.e., 2t = 2t+k = 2 <  for all t, k

24
The coefficient of autocorrelation (i.e. the ratio of covariance between Yt and Yt + k and variance of Yt) is also
used interchangeably.
25
For examples, the covariance between Yt and Yt-3 is the same the covariance between Yt-2 and Yt-5 (i.e. the
distance is 3 lags), but the covariance between Yt and Yt-3 is different from the covariance between Yt and Yt-5.
26
For more details, see Mills (2019: p.32-33).
27
Correlogram is a graph of the autocorrelations for various lags of a time series (Hanke and Wichern, 2014:
p.21). Autocorrelation is the correlation between a variable lagged one or more time periods and itself (Hanke and
Wichern, 2014: p.18).
20
 (c) covariance: Cov(Yt,Yt+k) = k = E[(Yt – )(Yt+k – )] (5.3)
i.e., Cov(Yt,Yt+k) = Cov(Yt+s,Yt+s+k) = k <  for all t, s, k
where k, covariance (or exactly autocovariance) at lag k, is the covariance between the
values of Yt and Yt-k, that is, between two Y values k periods apart. If k = 0, we obtain
0, which is simply the variance of Y (= 2); if k = 1, 1 is the covariance between two
adjacent values of Y. Note that we sometimes use Yt-k (where Yt is the current value of
Y) to replace Yt+k (where Yt is the the origin of Y). According to Asteriou and Hall
(2016: p.277), the quantities in above equations would remain the same whether the
observations for the series were, for example, from 2000 to 2010 or from 2010 to 2020.
Suppose we shift the origin of Y from Yt to Yt+m (say, from 1998Q3 to 2008Q3 for our
GDP data). Now, if Yt is to be stationary, the mean, variance, and autocovariance of
Yt+m must be the same as those of Yt. Translated into plain language, the above
properties means that a time series is (covariance) stationary, if its mean and variance
are constant and finite over time and if the covariance (or autocorrelation) depends only
on the time distance k between the two elements of the time series (i.e., the time between
the periods k) but not on the actual point in time t (Kocenda and Cerny, 2017: p.17-8;
Hill et al., 2018: p.427).
Gujarati and Porter (2009: p.741) state that such a time series will tend to return to its
mean (i.e., mean reversion)28 and fluctuations around its mean (measured by its
variance) will have a broadly constant value. In other words, a stationary process will
not drift far away from its mean value because of the finite variance. It also notes that
for a stationary process, the speed of mean reversion depends on the covariance; it is
quick if the covariances are small and slow when they are large (Gujarati and Porter,
2009: p.741). Shocks to a stationary series are necessarily temporary; over time, the
effects of the shocks will dissipate and the series will revert to its long-run mean level.
As such, long-term forecasts of a stationary series will converge to the unconditional
mean of the series (Asteriou and Hall, 2016: p.277). Another characteristic of stationary
variables is that they are weakly dependent, that is, their sample autocorrelations29 cut
off or tend to decline geometrically, dying out at long lags (Hill et al., 2018: p.566). In
other words, weak dependence implies that, as k →  (i.e., observations get further and
further apart in time), they become almost independent. For k large enough, the
autocorrelations become negligible (Hill et al., 2018: p.428).
If a time series is not stationary in the sense just defined, it is called a nonstationary
time series. In other words, a nonstationary time series will have a time-varying mean
or a time-varying variance or both (Gujarati and Porter, 2009: p.741). In addition,
another characteristic of nonstationary variables is that their sample autocorrelations
remain large at long lags, that is, they exhibit strong dependence (Hill et al., 2018:

28
In other words, if a time series is stationary, then any shock occurs in time t has a diminishing effect over time
and finally disappears in time t + k as k → . In contrast, the effect of a shock either remains present in the same
magnitude in all future dates (Kocenda and Cerny, 2017: p.17). Having a constant mean and fluctuations in the
series that tend to return to the mean are characteristics of stationary variables. This is called the property of mean
reversion (Hill et al., 2018: p.566).
29
Autocorrelation coefficients are introduced in the next section.

21
p.566). Unit root processes form a special subset of nonstationary processes, and they
play a particularly important role in time series analysis. Most economic time series are
not stationary and specific transformations are needed in order to achieve stationarity
(Gujarati and Porter, 2009: p.747; Kocenda and Cerny, 2017: p.18). Nonstationary
process and some useful transformations are described in the next section.
Why is stationarity important? There are at least two reasons. First, if a time series is
nonstationary, we can study its behaviour only for the time period under consideration
(Gujarati and Porter, 2009: p.741). Therefore, each set of time series data will be for a
particular episode. As a result, it is impossible to generalize it to other time periods.
Therefore, for the purpose of forecasting or policy analysis, such time series may have
little practical meaning. Second, if the series is nonstationary, all the typical results of
the classical ordinary least squares (OLS) regression are not valid. Regressions with
nonstationary series may be nonsense and are therefore called spurious (Gujarati and
Porter, 2009: p.748; Asteriou and Hall, 2016: p.277; Mills, 2019: p.237).
A special type of stationary stochastic process, namely, a purely random, or white noise
process, is frequently used in time series econometrics. A white noise is a time series
that does not contain any information that would help in estimation, except its variance
and higher moments (Kocenda and Cerny, 2017: p.17). According to Gujarati and Porter
(2009: p.741), we call a stochastic process purely random if it has zero mean, constant
variance 2, and is serially uncorrelated. Residuals from a correctly specified or a true
model that captures fully the data generating process are white noises (Kocenda and
Cerny, 2017: p.17). It is also noted that a stationary time series is not automatically
white noise. For white noise we need additional conditions that the mean and all
covariances30 are 0; i.e.,  = 0 and k = 0 for all lag k (Kocenda and Cerny, 2017: p.17).
In this notes, the white noise error process will be usually denoted as ut ~ iid(0,2); that
is, a series of identically and independently distributed random variable with 0 mean
and constant variance.

[[ 5.2 AR and MA Processes

In this section, we will investigate two typical types of the stationary process, namely
AutoRegressive (AR) and Moving Average (MA).

AR(1) process: A typical stationary process

The simplest, purely statistical time series model is the autoregressive of order one
model, or AR(1) model. This is defined as:
Xt = Xt-1 + ut, where ut ~ iid(0,2) (5.4)
where, for simplicity, we do not include a constant and || < 1, and ut is a Gaussian
(white noise) error term. According to Asteriou and Hall (2016: p.277), the assumption
behind the AR(1) model is that the time series behaviour of Xt is largely determined by
its own last value in the preceding period. So, what will happen in t is largely dependent

30
You will see in the next section that all autocorrelation coefficients (i.e. ACF) of a white noise are zero.

22
on what happened in t – 1. Alternatively, what will happen in t + 1 will be determined
by the behaviour of the series in the current time t.
For examples,
OilPt = 0.5OilPt-1 + ut
where OilPt-1 is change in the oil price in the last period and ut is any shock at the
current period.
If the constraint || < 1 is held, the series defined as equation (5.4) is stationary. Why?
If this is a case, the series Xt satisfies all three conditions for stationarity discussed in
Section 5.1 [Equations (5.1), (5.2), (5.3)].

Stationary properties of the AR(1) model: A proof

From equation (5.4), substituting repeatedly for lagged Xt, we have:

Xt = Xt-1 + ut, where ut ~ iid(0,2)
= [Xt-2 + ut-1] + ut
= 2Xt-2 + ut-1 + ut
= 2[Xt-3 + ut-2] + ut-1 + ut
= 3Xt-3 + 2ut-2 + ut-1 + ut
=…
= tX0 + 0ut-0 + 1ut-1 + 2ut-2 + … + t-2ut-(t-2) + t-1ut-(t-1)
= ut + ut-1 + 2ut-2 + … + t-2u2 + + t-1u1 (5.5)
Since || < 1, t will be close to zero for large t; and X0 is always zero because the
starting value of any time series is obviously zero. Thus we have that:
• Mean is constant (i.e., zero).
E[Xt] = E[ut + ut-1 + 2ut-2 + … + t-2u2 + t-1u1] = 0 (5.6)
In the same manner, we can prove that:
E(Xt-1) = E(Xt) = E(Xt+1) = 0 (5.7)
• Variance is constant.
Because of equation (5.5)
Xt = ut + ut-1 + 2ut-2 + … + t-2ut-(t-2) + t-1ut-(t-1)
= ut + Lut + 2L2ut + … + t-2Lt-2ut + t-1Lt-1ut
= (1 – L)-1ut as t → ∞ (5.8)31

31
According to infinite geometric progression:
S∞ = a + ar + ar2 + ar3 + ar4 + … + ar∞ where |r| < 1 (a)
rS∞ = ar + ar2 + ar3 + ar4 + ar5 + … + ar∞ (b)
(a) – (b) => S∞ – rS∞ = a => S∞ = a(1 – r)-1
=> S∞ = (1 – r)-1a (c)
23
So we have:
Xt-1 = ut-1 + ut-2 + 2ut-3 + … + t-2ut-(t-2)-1 + t-1ut-(t-1)-1
= ut-1 + Lut-1 + 2L2ut-1 + … + t-2Lt-2ut-1 + t-1Lt-1ut-1
= (1 – L)–1ut-1 as t → ∞ (5.9)30
And because
Var(ut) = Var(ut-1) = 2 (5.10)
Therefore,
Var(Xt) = Var(Xt-1) (5.11)
As a result, from equations (5.4), (5.10), and (5.11) we have:
Var(Xt) = Var(Xt-1) + Var(ut)
= 2Var(Xt-1) + Var(ut)
= 2Var(Xt) + Var(ut)
 (1 – 2)Var(Xt) = 2
 Var(Xt) = 2/(1 – 2) = constant if || < 1 (5.12)
• Covariance depends only on the time distance k.
From equation (5.4), we have
Xt+1 = Xt + ut+1 (5.13)
Xt+2 = Xt+1 + ut+2
= (Xt + ut+1) + ut+2
= 2Xt + ut+1 + ut+2 (5.14)
…
Xt+k = Xt+k-1 + ut+k
= [Xt+k-2 + ut+k-1] + ut+k
= 2Xt+k-2 + ut+k-1 + ut+k
= 2[Xt+k-3 + ut+k-2] + ut+k-1 + ut+k
= 3Xt+k-3 + 2ut+k-2 + ut+k-1 + ut+k
=…
= kXt+k-k + k-1ut+k-k-1 + … + 3ut+k-3 + 2ut+k-2 + ut+k-1 + ut+k
= kXt + ut+k + ut+k-1 + 2ut+k-2 + 3ut+k-3 + … + k-1ut-1
= kXt + vt (5.15)
Therefore, we have:
Cov(Xt, Xt+1) = E{[Xt – E(Xt)][Xt+1 – E(Xt+1)]}

In this case, Lut = ut-1, L2ut = ut-2, …, LTut = ut-T, L∞ut = ut-∞; L = r, and ut = a; and L is the lag operator.
24
 = E{[XtXt+1] – [YtE(Xt+1)] – [E(Xt)Xt+1] + [E(Xt)E(Xt+1)]}
= E[XtXt+1]
= E[Xt(Xt + ut+1)
= E(XtXt) + E(ut+1Xt)
= Var(Xt)
= 2/(1 – 2) (5.16)
Cov(Xt, Xt+2) = E{[Xt – E(Xt)][Xt+2 – E(Xt+2)]}
= E{[XtXt+2] – [YtE(Xt+2)] – [E(Xt)Xt+2] + [E(Xt)E(Xt+2)]}
= E[XtXt+2]
= E[Xt(2Xt + ut+1 + ut+2)
= E(2XtXt) + E(ut+1Xt) + E(ut+2Xt)
= 2Var(Xt)
= 22/(1 – 2) (5.17)
Cov(Xt, Xt+k) = E{[Xt – E(Xt)][Xt+k – E(Xt+k)]}
= E{[XtXt+k] – [YtE(Xt+k)] – [E(Xt)Xt+k] + [E(Xt)E(Xt+k)]}
= E[XtXt+k]
= E[Xt(kXt + vt)
= E(kXtXt) + E(vtXt)
= kVar(Xt)
= k2/(1 – 2) (5.18)
From (5.14) and (5.17), we have:
Therefore, the autocorrelation function32 of the AR(1) model will be given by:
k = corr(Xt, Xt+k) = Cov(Xt, Xt+k)/Var(Xt) = k (5.19)
where k = 1, 2, …, K or = –1, –2, …, –K. Since || < 1, k will decay exponentially as
k increases. We also conclude that the AR(1) model satisfies the weak dependence
property of a stationary series.

[ MA(1) process: A typical stationary process

The simplest moving average model is that of order one, or the MA(1) model, which is
defined as:

32
This is sometimes known as ‘internal correlation’ between a current value, X t, and previous, or lagged, values,
Xt-k, for k = 1, 2, … The lag-k (sample) autocorrelation is defined as:
∑𝑇𝑡=𝑘+1(𝑋𝑡 − 𝑋̅)(𝑋𝑡−𝑘 − 𝑋̅)
𝜌̂𝑘 =
𝑆2
where
𝑋̅ = 𝑇 −1 ∑𝑇𝑡=1 𝑋𝑡 and 𝑆 2 = 𝑇 −1 ∑𝑇𝑡=1(𝑋𝑡 − 𝑋̅)
The set of sample autocorrelations for various values of k is known as the sample autocorrelation function (ACF)
and plays a key role in time series analysis (Mills, 2019: p.3).
25
 Xt = ut + ut-1, where ut ~ iid(0,2) (5.20)
Therefore, the implication behind the MA(1) model is that Xt depends on the value of
the immediate past error, which is known at time t.
For examples,
OilPt = ut + 0.5ut-1
where OilPt is change in oil price and ut is typhoon at sea at the current time.
Lemonadet = ut - 0.5ut-1
where Lemonadet is change in lemonade quantity demanded and ut is change in
temperature at the current time.
The MA(1) model represents a stationary series because it satisfies all three conditions
for stationarity discussed in Section 5.1 [Equations (5.1), (5.2), (5.3)].

Stationary properties of the MA(1) model: A proof

From equation (5.20):

• Mean is constant.
E[Xt] = E[ut + ut-1] = E[ut] + E[ut-1] = 0 (5.21)
• Variance is constant.
Var(Xt) = Var(ut + ut-1) = Var(ut) + 2Var(ut-1)
= 2 + 22
= 2(1 + 2) (5.22)
Since both parameters  and  are constant, so Var(Xt) is indeed constant.
• Covariance depends only on the time distance k.
From equation (5.20), we have:
Xt+1 = ut+1 + ut (5.23)
Xt+2 = ut+2 + ut+1 (5.24)
Xt+k = ut+k + ut+k-1 (5.25)
So that we have:
Cov(Xt, Xt+1) = Cov(ut + ut-1, ut+1 + ut) = Cov(ut, ut) = 2 (5.26)
Cov(Xt+1, Xt+2) = Cov(ut+1 + ut, ut+2 + ut+1) = Cov(ut+1, ut+1) = 2 (5.27)
Cov(Xt, Xt+2) = Cov(ut + ut-1, ut+2 + ut+1) =0 (5.28)
Cov(Xt, Xt+k) = Cov(ut + ut-1, ut+k + ut+k-1) =0 k>1 (5.29)
In short, if the distance is 1, the covariance is different from 0; while the distance is
larger than 1, the covariance is always zero.
Therefore, the autocorrelation function of the MA(1) model will be given by:
θσ2 θ
= for k = 1
k = corr(Xt, Xt+k) = Cov(Xt, Xt+k)/Var(Xt) = { σ2 (1+θ2 ) (1+θ2 ) (5.30)
0 fork > 1

26
In other words, with the MA(1) model, the graph of the ACF is expected to have one
spike for k = 1, and then go down to zero immediately.
In sum, for AR(1) model, the autocorrelation function declines as k increases, while for
the MA(1) model, the autocorrelation function drops abruptly zero as k > 1.

AR(p) và MA(q) processes

A generalization of the AR(1) model is the AR(p) model; the number in parentheses
denotes the order of the autoregressive process and therefore the number of lagged
dependent variables in the model will have (Asteriou and Hall, 2016: p.279-80). For
example, the AR(2) model will be an autoregressive model of order two, and will have
the following form:
Xt = 1Xt-1 + 2Xt-2 + ut, where ut ~ iid(0,2) (5.31)
Similarly, the AR(p) model will be an autoregressive model of order p, and will have p
lagged terms, as in the following:
Xt = 1Xt-1 + 2Xt-2 + … + pXt-p + ut, where ut ~ iid(0,2) (5.32)
or, using the summation symbol:
p
Xt = ∑i=1 ∅i X t−i + ut (5.33)
Finally, using the lag operator L, we can write the AR(p) model as:
Xt(1 – 1L – 2L2 – … – pLp) = ut
(L)Xt = ut (5.34)
where (L)Xt is a polynomial function of Xt.
A necessary but not sufficient requirement for the AR(p) model to be stationary is that
the summation of the p autoregressive coefficients should be less than 1, that is:
∑pi=1 ∅i < 1 (5.35)
The general form of the moving average model is an MA(q) model of the form:
Xt = ut + 1ut-1 + 2ut-2 + … + qut-q, where ut ~ iid(0,2) (5.36)
which can be rewritten as:
q
Xt = ut + ∑j=1 i ut−j (5.37)
or, using the lag operator:
Xt = (1 + 1L + 2L2 + … = qLp)ut
= (L)ut (5.38)
Because any MA(q) process is, by definition, an average of q stationary white noise
processes, it follows that every moving average model is stationary, as long as q is finite
(Asteriou and Hall, 2016: p.283).

27
In principle, we can easily recognize whether a certain series follows AR(p) or MA(q)
by looking their corresponding partial autocorrelation functions (PACF) or
autocorrelation functions (ACF).

Autocorrelation functions

The ACF simply represents the pattern of autocorrelation coefficients with respect to
the corresponding lags. For example, AC1 is corr(Yt,Yt-1), AC2 corr(Yt,Yt-2), …, and
ACk corr(Yt,Yt-k).

γk covariance at lag k
ACk = k = = (5.39)33
γ0 variance
If k = 0, 0 = 1. Since both covariance and variance are measured in the same units of
measurements, k is a unitless, or pure, number (Gujarati, 2009: p.749). It lies between
– 1 and + 1, as any correlation coefficient does. If we plot k against k, the graph we
obtain is known as the population correlogram.
Since in practice we only have a realization (i.e., sample) of a stochastic process, we
can only compute the sample autocorrelation function (SACF), ρ̂𝑘 . To compute this, we
must first compute the sample covariance at lag k, γ̂𝑘 , and sample variance, γ̂0 , which
are defined as:
̅ )(Xt+k − X
∑(Xt − X ̅)
γ̂k = (5.40)
T−1

̅ )2
∑(Xt − X
γ̂0 = (5.41)
T−1
̅ is the sample mean.
where T is the sample size (i.e., number of realizations) and X
Therefore, the sample autocorrelation function at lag k is:
̂𝑘
γ
ρ̂𝑘 = ̂0
(5.42)
γ

which is simply the ratio of sample covariance (at lag k) to sample variance. A plot of
ρ̂𝑘 against k is known as the sample correlogram.

Partial autocorrelation functions

According to Gujarati and Porter (2009: p.778), the concept of partial autocorrelation is
analogous to the concept of partial regression coefficient. In the k-variable multiple
regression model, the kth regression coefficient Bk measures the rate of change in the
mean value of the regressand for a unit change in the kth regressor Xk, holding the
influence of all other regressors constant. In a simple language, the partial (sample)
regression coefficient of a three-regressor regression model can be defined as follows:
Y = b0 + b1X1 + b2X2 + b3X3 + e (5.43)

33
Due to stationary condition, Var(Xt-k) = Var(Xt) = Var(Xt+k), so the denominator of equation (5.39) is simply
Var(Yt) instead of SQRT[Var(Yt).Var(Yt+k)] (Kocenda and Cerny, 2017: p.42).
28
 X3 = c0 + c1X1 + c2X2 + v3 (5.44)
=> v3 = X3 – c0 – c1X1 – c2X2
(i.e., v3 is X4 after excluding the effects of X1 and X2)
Y = a0 + b3v3 + r (5.45)
Thus, b3 is the ‘direct’ effect of a unit change in X3 on the mean value of Y, net of any
effect that X1 and X2 may have on mean Y (see Appendix A for more details).
In the same token, the partial (sample) regression coefficient of a k-variable regression
model can be defined as follows:
Y = b0 + b1X1 + b2X2 + … + bkXk + e (5.46)
Xk = c0 + c1X1 + c2X2 + … + ck-1Xk-1 + vk (5.47)
=> vk = Xk – c0 – c1X1 – c2X2 – … – ck-1Xk-1
(i.e., vk is Xk after excluding the effects of other regressors on Xk)
Y = a0 + bkvk + r (5.48)
Thus, bk is the ‘direct’ effect of a unit change in Xk on the mean value of Y, net of any
effect that other regressors may have on mean Y.
In similar fashion, the PACF (k) measures correlation between observations that are k
periods apart after controlling for correlations at intermediate lags (i.e., lags less than
k). In other words, partial autocorrelation coefficient is the correlation between Yt and
Yt-k after removing the effect of the immediate Y’s. This means that in time series data,
a large proportion of correlation between Yt and Yt-k may be due to the correlations they
have with the intervening lags Yt-1, Yt-2, …, Yt-k+1. The partial autocorrelation
coefficient k removes the influence of these intervening variables. In a simple
language, the partial autocorrelation coefficient of a three-lag autoregressive model can
be defined as follows:
Yt = b0 + b1Yt-1 + b2Yt-2 + b3Yt-3 + et (5.49)
Yt-3 = c0 + c1Yt-1 + c2Yt-2 + v3 (5.50)
Yt = a0 + b3v3 + r (5.51)
Similarly, the partial (sample) autocorrelation coefficient of a k-lag autoregressive
model can be defined as follows:
Yt = b0 + b1Yt-1 + b2Yt-2 + … + bkYt-k + e (5.52)
Yt-k = c0 + c1Yt-1 + c2Yt-2 + … + ck-1Yt-k-1 + vk (5.53)
Yt = a0 + bkvk + r (5.54)
Thus, PAC1 is the regression coefficient of Yt-1 when regressing Yt on Yt-1; PAC2 is the
regression coefficient of Yt-2 when regressing Yt on Yt-1 and Yt-2; PAC3 is the regression
coefficient of Yt-3 when regressing Yt on Yt-1, Yt-2 and Yt-3, …, and PACp is the
regression coefficient of Yt-p when regressing Yt on Yt-1, Yt-2, …, and Yt-p. For example,
after executing the following do-file:
import excel using D:\Table13_6.xls, firstrow clear
29
 format date %td
tsset time
qui reg close [Link]
qui ereturn list
qui matrix b=e(b)
qui matrix list b
qui scalar PAC1=b[1,1]
qui reg close [Link] [Link]
qui ereturn list
qui matrix b=e(b)
qui matrix list b
qui scalar PAC2=b[1,2]
qui reg close [Link] [Link] [Link]
qui ereturn list
qui matrix b=e(b)
qui matrix list b
qui scalar PAC3=b[1,3]
qui reg [Link] [Link] [Link]
predict v, res
qui reg close v
qui ereturn list
qui matrix b=e(b)
qui matrix list b
qui scalar pac3=b[1,1]
display PAC1 PAC2 PAC3 pac3
You will see the results: PAC1 = 0.9906, PAC2 = 0.0537, PAC3 = pac3 = 0.0562.
The Stata’s corrgram command is usually used to identify both autocorrelation
function and partial autocorrelation function [corrgram variable, lag(k)].
corrgram close, lags(8)

30
To see both graph and statistics, we can use the following commands (but they are not
popularly used):
ac close, lags(12) generate(AC_close)
list AC_close in 1/12
pac close, lags(12) generate(PAC_close)
list PAC_close in 1/12

AR or MA process?

A rule of thumb

With a stationary series, the rule of thumb is as follows:

• AR(1) has PAC1  0, while the rest = 0; AR(2) has PAC1  0 and PAC2  0, while
the rest = 0; and so on.
• MA(1) has AC1  0, while the rest = 0; MA(2) has AC1  0 and AC2  0, while the
rest = 0; and so on.

A practical testing technique

How to know whether an individual AC or PAC coefficient is different from zero? The
most frequently used method is to look at the AC and PAC graphs (or the sample
correlogram). According to Hill et al. (2018: p.426), correlogram is a useful device for
assessing the significance of autocorrelation (and partial autocorrelation as well). An
autocorrelation (or a partial autocorrelation), rk, will be significantly different from zero
at a 5% significance level if SQRT(T)rk  1.96 or if SQRT(T)rk  –1.96. Alternatively,
we can say that rk will be significantly different from zero if rk  1.96/SQRT(T) or rk 
–1.96/SQRT(T). This is known as the 95% interval band. Where T is the sample size.
Experience tells that a certain coefficient that lies outside the 95% confidence band (i.e.,
the shaded area; but carefulness is needed if the series is nonstationary) appears to be
significantly different from zero at the 5% level, while it lies inside the 95% confidence
band may imply zero. More formally, we can use the t-statistic (in Section 8.2).
Approximate 95% confidence bands are computed using the fact that Z = (rk –
0)/SQRT(1/T) = SQRT(T)rk ~ N(0,1). Where rk is either ACk or PACk coefficients. In
Stata, we use the command gen z = sqrt(e(N))rk or gen z = sqrt(r(N))rk to generate the
boundary. If any of the numbers are less than –1.96 or greater than +1.96, it lies outside
of the approximate 95% confidence interval and it statistically significant at the 5%
level. Note that Stata’s ac and pac functions use a different method (Bartlett’s) and the
results may differ from those based on this simple approximation (Adkins and Hill,
2011: p.279).
import excel using D:\Table13_6.xls, firstrow clear
tsset time
ac close, lags(12) gen(ac_close)
qui sum close
qui return list
31
 gen z = sqrt(r(N))*ac_close
list ac_close z in 1/10
qui pac close, lags(12) gen(pac_close)
qui sum close
gen Z = sqrt(r(N))*pac_close
list pac_close Z in 1/10

5.3 Examples of MA and AR Processes

MA(1) and AR(1)

Using Stata do-file as the following commands:

clear
set obs 400
gen timevar = _n
set seed 12345
drawnorm e, n(400) means(0) sds(1)
tsset34 timevar
gen AR1 = 0
replace AR1 = 0.7*L.AR1 + e if _n > 1
gen MA1 = 0
replace MA1 = e+0.7*L.e if _n > 1
ac MA1, lags(15)
pac MA1, lags(15)

34
Again, we must use the command ‘tsset’ to declare that the dataset is a time series data. If we ignore this
command, other commands such as tlline, pac, dfuller, … cannot be performed correctly.
32
 We have the following graphs:

0.40
0.40
0.30

Partial autocorrelations of MA1

0.20
Autocorrelations of MA1

0.20

0.00
0.10

-0.20
0.00
-0.10

-0.40
0 5 10 15 0 5 10 15
Lag Lag
Bartlett's formula for MA(q) 95% confidence bands 95% Confidence bands [se = 1/sqrt(n)]

Figure 5.1: ACF of MA(1) process. Figure 5.2: PACF of MA(1) process.

We see that only AC1 (i.e., 1) of MA(1) process is significantly different from zero at
the 5% level (i.e., outside the shaded area) and other coefficients of MA(1) process (i.e.,
2, 3, …, 15) are zero (inside the shaded area).
ac AR1, lags(15)
pac AR1, lags(15)
We have the following graphs:
0.80

0.80
0.60

0.60
Partial autocorrelations of AR1
Autocorrelations of AR1

0.40

0.40
0.20

0.20
0.00

0.00
-0.20

-0.20

0 5 10 15 0 5 10 15
Lag Lag
Bartlett's formula for MA(q) 95% confidence bands 95% Confidence bands [se = 1/sqrt(n)]

Figure 5.3: ACF of AR(1) process. Figure 5.4: PACF of AR(1) process.

We see that ACk (i.e., k = k) of AR(1) process is declining to zero as k increases. In
this case, only PAC1 (i.e., 1) is significantly different from zero at the 5% level (i.e.,

33
 outside the shaded area) and other coefficients of AR(1) process (i.e., 2, 3, 4,…, 15)
are zero (i.e., inside the shaded area).

MA(2) and AR(2)

Using Stata do-file as the following commands:

clear
set obs 400
gen t = _n
tsset t
gen eps = invnorm(uniform())
scalar theta0= 0
scalar theta1= 0.6
scalar theta2= 0.3
gen double AR2 = 0
qui replace AR2 in 3/l = theta0 + theta1*L.AR2 + theta2*L2.AR2 + eps
gen double MA2 =0
qui replace MA2 in 3/l = eps + theta1*[Link] + theta2*[Link]
ac MA2, lags(15)
pac MA2, lags(15)
We have the following graphs:
0.60

0.60
Partial autocorrelations of MA2
0.40

0.40
Autocorrelations of MA2

0.20

0.20
0.00

0.00
-0.20

-0.20

0 5 10 15 0 5 10 15
Lag Lag
Bartlett's formula for MA(q) 95% confidence bands 95% Confidence bands [se = 1/sqrt(n)]

Figure 5.5: ACF of MA(2) process. Figure 5.6: PACF of MA(2) process.

We see that only AC1 and AC2 (i.e., 1 and 2) of MA(2) process are statistically
different from zero (i.e., outside the shaded area) and other coefficients of MA(2)
process (i.e., 3, 4, …, 15) are zero (i.e., inside the shaded area).
ac AR2, lags(15)

34
 pac AR2, lags(15)
We have the following graphs:
0.80

0.80
0.60

Partial autocorrelations of AR2

0.60
Autocorrelations of AR2

0.40

0.40
0.20

0.20
0.00

0.00
-0.20

0 5 10 15 -0.20 0 5 10 15
Lag Lag
Bartlett's formula for MA(q) 95% confidence bands 95% Confidence bands [se = 1/sqrt(n)]

Figure 5.7: ACF of AR(2) process. Figure 5.8: PACF of AR(2) process.

We see that ACk (i.e., k) of AR(2) process is declining to zero as k increases. In this
case, only PAC1 and PAC2 (i.e., 1 and 2) are statistically different from zero (i.e.,
outside the shaded area) and other coefficients of AR(2) process (i.e., 3, 4, …, 15) are
zero (i.e., inside the shaded area).
By the same token, you can generate other series with higher lag orders such as MA(3)
and AR(3), MA(4) and AR(4), and so on. It is important to note that the AR(p) model
is stationary if the condition in equation (30) is satisfied. It is also important to know
that actual economic series rarely exhibit the exact patterns as theoretically shown.

5.4 Invertibility: Converting an AR(1) to an MA() and Vice Versa

According Asteriou and Hall (2016: p.283), a property often discussed in connection
with the moving average processes is that of invertibility. A time series Xt is invertible
if it can be represented by a finite-order moving average or convergent autoregressive
process. Invertibility is important because the use of the ACF and PACF for
identification assumes implicitly that the Xt sequence can be approximated well by an
autoregressive model.

AR(1) => MA()

If || < 1, the AR(1) process can be converted into an infinite order MA process with
geometrically declining weights. This is simply proved as follows:
Xt = Xt-1 + ut, where ut ~ iid(0,2) (5.4)
= [Xt-2 + ut-1] + ut
= 2Xt-2 + ut-1 + ut
35
 =…
= Xt- + 0ut-0 + 1ut-1 + 2ut-2 + … + ut-
= 0ut-0 + 1ut-1 + 2ut-2 + … + ut-
= ut + ut-1 + 2ut-2 + … + ut-
= ut + 1ut-1 + 2ut-2 + … + ut- = MA() (5.55)35

MA(1) => AR()

If || < 1, the MA(1) process can be converted into an infinite order AR process with
geometrically declining weights. This is simply proved as follows:
Xt = ut – ut-1, where ut ~ iid(0,2) (5.56)36
Using the lag operator [i.e., Lut = ut-1, L2ut = ut-2, L3ut = ut-3, …], equation (5.56) can be
rewritten as:
Xt = (1 – L)ut
Xt
ut = (5.57)
(1− θL)

If || < 1, then the left-hand side of equation (5.57) can be considered as the sum of an
infinite geometric progression:
ut = Xt(1 + L + 2L2 + 3L3 + …)
ut = Xt + LXt + 2L2Xt + 3L3Xt + …
Xt = – LXt – 2L2Xt – 3L3Xt – … + ut
Xt = – Xt-1 – 2Xt-2 – 3Xt-3 – … + ut
Xt = 1Xt-1 + 2Xt-2 + 3Xt-3 + … + ut (5.58)
To understand equation (5.58), let us rewrite the MA(1) model as defined in equation
(5.56) as below:
ut = Xt + ut-1 (5.59)
Lagging the relationship in equation (5.59) one period, we have:
ut-1 = Xt-1 + ut-2 (1.60)
Substituting this into the original expression [i.e., equation (5.56)], we have:
Xt = ut – [Xt-1 + ut-2] = ut – Xt-1 – 2ut-2 (5.61)
Lagging the relationship in equation (5.59) one period, we have:
ut-2 = Xt-2 + ut-3 (5.62)

35
Another way:
Xt = LXt + ut => (1 – L)Xt = ut => Xt = (1 – L)-1ut = (1 + L + 2L2 + 3L3 + …)ut = ut + ut-1 + 2ut-2 + 3ut-3
+ … if || < 1 (Mills, 2019: p.35).
36
We can apply for the case in equation (5.20): Xt = ut + ut-1 (see Asteriou and Hall, 2016: p.283-4).
36
Substituting equation (5.62) into the expression in equation (5.61), we have:
Xt = ut – Xt-1 – 2[Xt-2 + ut-3] = ut – Xt-1 – 2Xt-2 – 3ut-3 (5.63)
If we go on this procedure [i.e., lagging and substituting] for an infinite number of times,
we finally get the expression in equation (5.58).

[[ 5.5 ARMA Process

The model

ARMA(p,q) models are combinations of AR(p) and MA(q) processes. The general form
of the ARMA model is defined as:
Xt = 1Xt-1 + 2Xt-2 + … + pXt-p + ut
+ 1ut-1 + 2ut-2 + … + qut-q, where ut ~ iid(0,2) (5.64)
which can be rewritten, using the summation, as:
p q
X t = ∑i=1 θi X t−i + ut + ∑j=1 ϕj ut−j (5.65)
or, using the lag operator:
Yt (1 − ϕ1 L − ϕ2 L2 − ⋯ − ϕp Lp ) = (1 + θ1 L + θ2 L2 + ⋯ + θq Lq )ut
Φ(𝐿 ) X t = Θ (𝐿 )u t (5.66)
According to Asteriou and Hall (2016: p.285), in the ARMA(p,q) models the condition
for stationarity deals only with the AR(p) part of the specification. Similarly, the
property of invertibility for the ARMA(p,q) models will relate only with the MA(q) part
of the specification. It is also important to note that the ARMA(p,q) models are specified
for the stationary series. Therefore, it the series is not stationary, the series needs to be
transformed to make it stationary (i.e., often take logarithms and then first differences)
before applying the ARMA(p,q) models. This refers to the ARIMA(p,d,q) models. The
ARIMA models and Stata’s commands are discussed in Appendix B.

Identification

Once stationarity has been achieved, it is necessary to identify the p and q orders of the
ARMA model. For a pure MA(q) process, the ACF will tend to show estimates that are
significantly different from zero up to lag q and then die down immediately after qth lag.
The PACF for MA(q) will tend to die down quickly, either by an exponential decay or
by a damped sine wave. In contrast, the pure AR(p) process will have an ACF that will
tend to die down quickly, either by an exponential decay or by a damped sine wave,
while the PACF will tend to shows spikes (i.e., significant partial autocorrelation) for
lags up to p and the will die down immediately (Asteriou and Hall, 2016: p. 287).
If neither the ACF nor the PACF show a definite cut-off, a mixed process is suggested.
In this case it is difficult, but not possible, to identify the AR and MA orders. We should
think of the ACF and PACF of pure AR and MA processes as being superimposed onto
one another. For example, if both ACF and PACF shows signs of slow slow exponential
decay, an ARMA(1,1) process may be identified. Similarly, if the ACF shows three
37
significant spikes at lags one, two and three and then an exponential decay, and the
PACF spikes at the first lag and then shows an exponential decay, an ARMA(1,3). Table
5.1 and Figure 5.9 report some possible combinations of ACF and PACF forms that
allow us the detection of the order of ARMA processes.
Table 5.1: ACF and PACF patterns for ARMA(p,q) models.
Model ACF PACF
MA(1) Single significant spike at lag 1 Damped sine wave or exponential decay
AR(1) Damped sine wave or exp. decay Single significant spike at lag 1
MA(2) Significant spikes up to lag 2 Damped sine wave or exp. decay
AR(2) Damped sine wave or exp. decay Significant spikes up to lag 2
MA(q) Significant spikes up to lag q Damped sine wave or exp. decay
AR(p) Damped sine wave or exp. decay Significant spikes up to lag p
ARMA(1,1) Decay (exp. or sine wave) Decay (exp. or sine wave) beginning at
beginning at lag 1 lag 1
ARMA(2,2) Decay (exp. or sine wave) Decay (exp. or sine wave) beginning at
beginning at lag 2 lag 2
ARMA(p,q) Decay (exp. or sine wave) Decay (exp. or sine wave) beginning at
beginning at lag q lag p

Figure 5.9: ACF and PACF patterns for ARMA(p,q) models.

Model ACF PACF

(a)
AR(1)

(a)

(b) (b)

38
AR(2)

(c) (c)

(d) (d)

(e) (e)
MA(1)

(f) (f)

39
 (g) (g)
MA(2)

(h) (h)

(i) (i)
ARMA(1,1)

(j) (j)

(k) (k)

40
 (l) (l)
ARMA(2,2)

(m) (m)

(n) (n)

Source: Hanke and Wichern (2014: p.357-9).

41
 6.
NONSTATIONARY STOCHASTIC PROCESSES

6.1 Trends and Random Walk

According to Stock and Watson (2015: p.523), time series variables can fail to be
stationary in various ways, but two are especially relevant for regression analysis of
economic time series: (1) the series can have persistent, long-run movements, that is,
the series can have trends; and, (2) the population regression can be unstable over time,
that is, the population regression can have breaks.
For the purpose of this series of lectures, we mainly focus on the first type of non-
stationarity. Nonstationary variables that wander up and down, trending in one direction
and then the other, are said to posses a stochastic trend. A trend is a persistent long-term
movement of a variable over time37. A time series variable fluctuates around its trend38.
There are two types of trends often observed in time series data: deterministic or/and
stochastic (Hill et al., 2018: p.567). It may be noted that a time series may contain both
a stochastic and a deterministic trend (Gujarati, 2009: p.575).

Trend stationary variables

Variables that are stationary after “subtracting” a deterministic trend39 are called trend
stationary (Hill et al., 2018: p.567). The simplest model for a deterministic trend for a
variable Y is the linear trend model:
Yt = A + Bt + ut, where ut ~ iid(0,2) (6.1)
where t = 1, 2, …, T; t is known as a trend or time variable. If we focus on the trend and
assume any change in the error is zero (i.e., ut = ut – ut-1 = 0), then the coefficient 
gives a change in Y from a one time period to the next:
Yt – Yt-1 = (A + Bt) – |A + B(t–1)| + ut = B (6.2)
The trend A + Bt is called a deterministic trend because it does not contain a stochastic
(random) component. The variable Yt is trend stationary if its fluctuations around this
trend are stationary. Since these fluctuations are given by changes in the error term:
ut = Yt – (A + Bt) (6.3)
is called a “detrended variable”. Yt is trend stationary if ut is stationary. According to
Hill et al. (2018: p.568), when Yt is trend stationary, we can use OLS to find estimates

37
The trend is the long-term component that represents the growth or decline in the time series over an extended
period of time (Hanke and Wichern, 2014: p.17).
38
When observations exhibit rises and falls that are not of a fixed period, a cyclical pattern exists. The cyclical
component is the wavelike fluctuation around the trend that is usually affected by general economic conditions
(Hanke and Wichern, 2014: p.17).
39
Mills (2019: p.6) calls this by the terms ‘linear trend’ and ‘non-linear trend’.
42
a and b from equation (6.1) and then convert the trend stationary variable Yt to a
stationary variable et by removing the trend:
et = Yt – (a + bt) (6.4)
It is useful to note that if we are considering a regression or an autoregressive distributed
lag (ARDL)40 model involving two trend stationary variables, say Yt and Xt, we must
use the detrended variables instead of using the original ones (Hill et al., 2018: p.568).
Such models are discussed in Section 9. It is also useful to emphasize that the trend
stationary variables are often referred to as “stationary”, although they are not actually
stationary because their means change over time (i.e., E(Yt) = Bt). In addition to the
linear trend stationary model in (62), others are possible. For examples, we might have
to use either a quadratic trend (i.e., Yt = A + Bt + Ct2 + ut) or a cubic trend (i.e., Yt = A
+ Bt + Ct2 + Dt3 + ut) depending on the data pattern. In some cases, we might use the
exponential trend (i.e., Yt = e(A + Bt + ut) => ln(Yt) = A + Bt + ut). The variable ln(Yt) is
trend stationary if ut is stationary (Hill et al., 2018: p.569).
The Stata’s commands for creating a detrended variable:
clear
set obs 100
gen quarter = tq(1995q1) + _n-1
format quarter %tq
set seed 12345
drawnorm u, n(100) means(0) sds(1)
gen time = _n
gen Y = time + u
reg Y time
predict e, residuals
tsset quarter
tsline Y
tsline e
3
6

2
4

1
Residuals
2
Y

0
0

-1
-2
-4

-2

1995q1 2000q1 2005q1 2010q1 2015q1 2020q1 1995q1 2000q1 2005q1 2010q1 2015q1 2020q1
quarter quarter

Figure 6.1: Linear trend stationary variable. Figure 6.2: Detrended variable.

40
See Appendix C.
43
Random walk models and difference stationary variables

In contrast, a stochastic trend is random and varies over time. According to Stock and
Watson (2015: p.552), it is more appropriate to model economic time series as having
stochastic rather than deterministic trends. Therefore, our treatment of trends in
economic time series data focuses mainly on stochastic rather than deterministic trends,
and when we refer to “trends” in time series data, we mean stochastic trends unless we
explicitly say otherwise.
The simplest model of a variable with a stochastic trend is the random walk41. There are
two types of random walks: (1) random walk without drift (i.e., no constant or intercept
term) and (2) random walk with drift (i.e., a constant term is present).

Random walk without drift

The random walk without drift is defined as follow. Suppose ut is a white noise error
term with mean 0 and variance 2. The Yt is said to be a random walk without drift if:
Yt = Yt-1 + ut, where ut ~ iid(0,2) (6.5)
This equation is just a special case of equation (5.4), except  is now equal 142. In
statistical language, ‘ = 1’ is called ‘unit root’. The basic idea of a random walk is that
the value of the series tomorrow (Yt+1) is its value today (Yt), plus an unpredictable
change (ut+1).
We can understand the “wandering” behaviour of random walk models by doing some
recursive substitution. From equation (6.5), we can write:
Y1 = Y 0 + u 1
Y2 = Y 1 + u 2 = Y 0 + u 1 + u 2
Y3 = Y 2 + u 3 = Y 0 + u 1 + u 2 + u 3
…
Yt = Yt-1 + ut = Y0 + u1 + … + ut (6.6)
Therefore, the random walk model contains an initial value Y0 (often set to zero because
it is so far in the past that its contribution to Yt is negligible) plus a component that is
the sum of the past stochastic terms ∑Tt=1 ut . This latter component is thus called the
stochastic trend. This term arises because a stochastic component ut is added for each
time t, and because it causes the time series to trend in unpredictable directions. If the
variable Yt is subjected to a sequence of positive shocks, ut > 0, followed by a sequence
of negative shocks, ut < 0, it will have the appearance of wandering upward, then
downward (Hill et al., 2018: p.573).
What are the properties of the random walk without drift? If the process started at some
time 0 with a value Y0 = 0, equation (6.6) can be rewritten as:
Yt = Y0 + ∑Tt=1 ut (6.7)

41
The time series is called a random walk because it appears to wander slowly upward or downward with no real
pattern (Hill et al., 2018: p.573).
42
Again, remember that if || < 1, Yt is an AR(1) process, and Yt is stationary.
44
 therefore,
E(Yt ) = E(Y0 + ∑Tt=1 ut ) = Y0 = 0 (6.8)
In like fashion, it can be shown that
Var(Yt) = E(Y0 + ∑ ut – Y0)2 = E(∑ ut )2 = t2 (6.9)
Therefore, the random walk without drift has a mean equal to its initial or starting value
(which is constant), and a variance that increases over time, eventually becoming
infinite. Although the mean is constant, the increasing variance implies that the series
may not return to its mean, and so sample means taken for different periods are not the
same. In other words, the variance of Yt depends on t, then its distribution depends on
t; that is, it is nonstationary. Interestingly, if we re-write equation (6.5) as:
(Yt – Yt-1) = ∆Yt = ut (6.10)
where ∆Yt is the first difference of Yt. It is easy to show that, while Yt is nonstationary,
its first difference is stationary (why?)43. And this is very significant when we work
with time series data. For examples, unit root tests and most modelling techniques of
time series data are based on the first differences. In its terminology, random walks are
widely known as the difference stationary variables. Differencing is an important data
transformation method in time series analysis.
Using Stata do-file as the following commands:
clear gen RW = 0
set obs 500 replace RW = [Link] + e if _n > 1
gen time = _n label variable RW "Random walk ||
set seed 12345 without drift"
drawnorm e, n(500) || tsline RW
means(0) sds(1) tsline [Link]
tsset time
We have the following graphs:
40

4
Random walk without drift, D
30

2
Random walk without drift

20

0
10

-2
-4
0

0 100 200 300 400 500 0 100 200 300 400 500
time time

Figure 6.3: Random walk without drift. Figure 6.4: ∆ random walk without drift.

43
Because ut is known as the white noise, which is assumed to follow a normal distribution with zero mean and
constant variance.
45
Random walk with drift

The random walk with drift can be defined as follow:

Yt =  + Yt-1 + ut, where ut ~ iid(0,2) (6.11)
where the constant term  is known as the drift parameter. Equation (6.11) shows that
each realization of the random variable Yt contains an intercept (the drift component )
plus last period’s value Yt-1 plus the error ut. According to Hill et al. (2018: p.573), the
time series data of this type appear to be “wandering” as well as “trending” upward or
downward. In general, random walk with drift models show definite trends either
upward (when the drift  is positive) or downward (when the drift  is negative).
By applying recursive substitution, we can show that the random walk with drift also
violates conditions of stationarity. Equation (6.11) can be rewritten as:
Y1 =  + Y0 + u 1
Y2 =  + Y1 + u 2 =  +  + Y0 + u 1 + u 2
Y3 =  + Y 2 + u 3 =  +  +  + Y 0 + u 1 + u 2 + u 3
…
YT = YT-1 + uT =  +  + … +  + Y0 + u1 + … + uT
Yt = t + Y0 + ∑Tt=1 ut (6.12)
Equation (6.12) shows that the value of Y at time t is made up of an initial value Y0, the
stochastic trend component ∑Tt=1 ut , and now a deterministic trend component t.
According to Hill et al. (2018: p.573-4), the component t is called a deterministic trend
because a fixed value  is added for each time t. The mean and variance of Yt are:
E(Yt) = Y0 + t (6.13)
Var(Yt) = t2 (6.14)
In other words, both mean and variance of Yt depends on t, its distribution depends on
t, that is, it is nonstationary.
We can extend the random walk model even further by adding a time trend:
Yt =  + t + Yt-1 + ut, where ut ~ iid(0,2) (6.15)
Note that the addition of a time trend variable t strengthens the trend behaviour (Hill et
al., 2018: p.574). We can see the amplification using the same algebraic manipulation
as before:
Y1 =  +  + Y0 + u 1
Y2 =  + 2 + Y1 + u2 =  + 2 + [ +  + Y0 + u1] + u2
= 2 + 3 + Y0 + u1 + u2
Y3 =  + 3 + Y2 + u3 =  + 3 + [2 + 3 + Y0 + u1 + u2] + u3
= 3 + 6 + Y0 + u1 + u2 + u3
…
t(t+1)
Yt = t + (
2
)  + Y0 + ∑Tt=1 ut (6.16)

46
Therefore, the additional term has the effect of strengthening the trend behavior. Similar
to the random walk with drift, both mean and variance of Yt depends on t, its distribution
depends on t, that is, it is nonstationary. In such a case, we have to remove the time
trend of the first difference in order to transform to the stationary variable.
Using Stata do-file as the following commands:
clear
set obs 500
gen time = _n
set seed 12345
drawnorm e, n(500) means(0) sds(1)
tsset time
gen RW = 0
replace RW = 0.2 + [Link] + e if _n > 1
label variable RW "Random walk with a positive drift"
tsline RW
tsline [Link]
We have the following graphs:
150

4
Random walk with a positive drift, D

2
100

0
50

-2
-4
0

0 100 200 300 400 500 0 100 200 300 400 500
time time

Figure 6.5: Random walk with a positive drift. Figure 6.6: ∆ random walk with a positive drift.

Using Stata do-file as the following commands:

clear
set obs 500
gen time = _n
set seed 12345
drawnorm e, n(500) means(0) sds(1)
tsset time
gen RW = 0
replace RW = -0.2 + [Link] + e if _n > 1
label variable RW "Random walk with a negative drift"
tsline RW
tsline [Link]
47
We have the following graphs:

4
0

Random walk with a negative drift, D

2
-20

0
-40

-2
-60
-80

-4
0 100 200 300 400 500 0 100 200 300 400 500
time time

Figure 6.7: A random walk with a negative drift. Figure 6.8: ∆ random walk with a negative drift.

Using Stata do-file as the following commands:

clear
set obs 500
gen time = _n
set seed 12345
drawnorm e, n(500) means(0) sds(1)
tsset time
gen RW = 0
replace RW = 0.2 + 0.01*time + [Link] + e if _n > 1
label variable RW "Random walk with a time trend"
tsline RW
tsline [Link]
reg [Link] time
predict detrend, residuals
label variable detrend “First difference after detrending”
tsline detrend
1500

8
Random walk with a time trend, D
6
Random walk with a time trend
1000

2 4
500

0 -2
0

0 100 200 300 400 500 0 100 200 300 400 500
time time

Figure 6.9: Random walk with a time trend. Figure 6.10: ∆ random walk with a time trend.

48
 4
First difference after detrending
-2 0-4 2

0 100 200 300 400 500

time

Figure 6.11: ∆ random walk with a time trend after detrending.

In a nutshell, a random walk is a nonstationary process, where either its mean or its
variance or both increases over time. Because stochastic trends are the most prevalent
source of non-stationarity, and they introduce special problems, when we refer to
nonstationary variables, we will generally mean variables that are neither stationary nor
trend stationary. In other words, we often imply difference stationary variables or
nonstationary variables with stochastic trends (Hill et al., 2018: p.576).

Trend versus difference stationary: An example

Using Stata do-file as the following commands:

use "D:\Table13_1.dta", clear
tsset time
label variable lnex "Log of the dollar/euro daily exchange rate"
tsline lnex
reg lnex time
predict e, residual
label variable e "detrended variable of log(ex)"
tsline e
tsline [Link]
For example, the lnex [the logarithm44 of the dollar/euro daily exchange rate, i.e., lnex
= log(ex), see data on [Link], Gujarati (2011)45] is a nonstationary series (Figure
6.12), and its detrended series (i.e., residuals from the OLS regression of lnex on time

44
Taking a natural logarithm is applied when the data perform exponential growth, which is a common case in
economics. Prior to differencing and detrending, the most common transformation is to take a natural logarithm
of the data in order to deal with a sort of non-linearity or to reduce an exponential trend into a linear one. In
addition, if we are interested in growth rates, it is natural to apply logarithmic differencing (Kocenda and Cerny,
2017: p.19).
45
This dataset can be downloaded here: [Link]
Example/student-zone/.
49
variable: et = lnext – a – bt) is still nonstationary (Figure 6.13). This points out that the
variable lnex is not a trend stationary variable. Figures 6.14 shows that the logarithm of
the dollar/euro daily exchange rate is characterized as a difference stationary process
because its level is not stationary, whereas its first difference is stationary. Formally,
we will use statistical tests for non-stationarity, called unit root tests such as Dickey-
Fuller and Phillips-Perron. These tests will be discussed in Section 8.
.6

.3 .2
detrended variable of log(ex)
.4

.1
.2

0
0

-.1
-.2

-.2

0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
Time Time

Figure 6.12: Log of the $/€ daily exchange rate. Figure 6.13: Detrended variable of log(EX).
.04
Log of the dollar/euro daily exchange rate, D
-.02 0 .02

0 500 1000 1500 2000 2500

Time

Figure 6.14: First difference of log(EX).

50
6.2 Unit Root Stochastic Process
The simplest model

According to Gujarati and Porter (2009: p.744), the random walk model is an example
of what is known in the literature as a unit root process.
Let us write the random walk model (6.15) as:
Yt = Yt-1 + ut, where ut ~ iid(0,2) and (–1    1) (6.17)
This model resembles the Markov first-order autoregressive model [AR(1)], as
mentioned in equation (5.4). If  = 1, equation (6.17) becomes a random walk without
drift. If  is in fact 1, we face what is known as the unit root problem, that is, a situation
of non-stationarity. The name unit root is due to the fact that  = 1. Technically, if  =
1, equation (6.17) becomes equation (6.15), then Yt – Yt-1 = ut. Now using the lag
operator L so that Lyt = Yt-1, L2Yt = Yt-2, L3Yt = Yt-3, and so on, we can write equation
(6.15) as (1 – L)Yt = ut. If we set (1 – L) = 0, we obtain, L = 1, hence the name ‘unit
root’. Thus, the terms non-stationarity, random walk, and unit root can be treated as
synonymous.
If, however, || < 1 (usually  < 1), that is if the absolute value of  is less than one, then
it can be shown that the time series Yt is stationary. In other words, if  is less than 1,
the equation (6.17) is an AR(1) process, which is proved as a stationary process with
constant mean, constant variance, and time-invariant covariance (see Section 5.2). In
sum, the series with  = 1 is a random walk (or nonstationary series), while the series
with || < 1 is an AR(1) process (or stationary series).

Illustrative examples

Consider the AR(1) model as presented in equation (6.17). Generally, we can have three
possible cases:
Case 1: || < 1 and therefore the series Yt is stationary. Graphs of stationary series for
 = 0.67 and  = 0.98 are presented in Figures 6.15 and 6.16.
Case 2:  = 1 where in this case the series contains a unit root and is non-stationary.
Graph of stationary series for  = 1 is presented in Figure 6.17.
Case 3: || > 1 where in this case the series explodes. A graph of an explosive series
for  = 1.26 is presented in Figure 6.18.
In order to reproduce the graphs and the series which are stationary, exploding and
nonstationary, we type the following commands in Stata:
clear
set obs 500
gen timevar = _n
set seed 12345
tsset timevar
drawnorm e1 e2 e3 e4, n(500) means(0 0 0 0) sds(1 1 1 1)
gen AR_67 = 0
51
 replace AR_67 = 0.67*L.AR_67 + e1 if _n > 1
label variable AR_67 "AR(1) with theta = 0.67"
gen AR_97 = 0
replace AR_97 = 0.97*L.AR_97 + e2 if _n > 1
label variable AR_97 "AR(1) with theta = 0.97"
gen RW = 0
replace RW = [Link] + e3 if _n > 1
label variable RW "AR(1) with theta = 1"
gen EP = 0
replace EP = 1.1*[Link] + e4 if _n > 1
label variable EP "AR(1) with theta = 1.1"
tsline AR_67
tsline AR_97
tsline RW
tsline EP

10
6

5
4

AR(1) with theta = 0.97

AR(1) with theta = 0.67

0
2

-5
0

-10
-2

-15
-4

0 100 200 300 400 500 0 100 200 300 400 500
timevar timevar

Figure 6.15: AR(1) with  = 0.67 (stationary). Figure 6.16: AR(1) with  = 0.97 (stationary).
4.00e+19
60

3.00e+19
AR(1) with theta = 1.1
AR(1) with theta = 1

40

2.00e+19
20

1.00e+19
0

0 100 200 300 400 500 0 100 200 300 400 500
timevar timevar

Figure 6.17: AR(1) with  = 1 (random walk). Figure 6.18: AR(1) with  = 1.1 (explosive).

52
 7.
SPURIOUS REGRESSION

7.1 Understanding the Concept

A consequence of performing the OLS regression involving nonstationary variables is
that OLS estimates no longer have approximate normal distributions in even large T.
That means interval estimates and hypothesis tests will no longer be valid. Precision of
estimation may not be what it seems to be and conclusions about relationships between
variables could be incorrect. One particular problem is that two completely independent
random walks can appear to have a strong linear relationship when none exists. This
phenomenon is known as spurious regression46. In such a case the OLS estimators do
not have their usual properties, and t-statistics are not reliable (Hill et al., 2018: p.574-
5). Specifically, the regression results from nonstationary variables often have very high
value of R2 (sometimes even higher than 0.95) and very high values of t-statistics
(sometimes even higher than 4), although these variables might not have any real
interrelationships (Asteriou and Hall, 2016: p.351). Since many macroeconomic time
series are nonstationary, it is particularly important to take care when estimating
regressions with macroeconomic variables (Hill et al., 2018: p.575).
What is an essence of a spurious regression? The intuition behind this is quite simple.
Over time, we expect any nonstationary series to wander around, so over any reasonably
long sample the series either drift up or down. If we then consider two completely
unrelated series which are both nonstationary, we would expect that either they will
both go up or down together, or one will go up while the other goes down. If we run a
OLS regression of one series on another, we would then find either a significant positive
relationship if they are going in the same direction or a significant negative one if they
are going in opposite directions even though they are really unrelated.
Granger and Newbold (1974) constructed a Monte Carlo analysis47 generating a large
number of Yt and Xt series containing unit roots following the formulas:
Yt = Yt-1 + uYt (7.1)
Xt = Xt-1 + uXt (7.2)

46
Granger and Newbold (1974) argued that regressions using integrated time series were likely to be spurious, in
that they would typically produce an apparently significant relationship even when the variables were unrelated
to each other (Mills, 2019: p.237). In addition, the spurious regression problem is not resolved by attempting to
remove deterministic trends from the data (Mills, 2019: p.241). According to Greene (2018: p.1026), in a seminal
paper, Granger and Newbold (1974) argued that researchers had not paid sufficient attention to the warning of
very high autocorrelation in the residuals from conventional regression models. The standard significance tests
such as t and F tests from regressions involving the levels of macroeconomic data are usually misleading, that is,
tend to reject the hypothesis of no relationship when there might be none. Granger and Newbold (1974) also
suggested that researchers should use a critical t value of 11.2 rather than the standard normal value of 1.96 to
assess the significance of a coefficient estimate.
47
This was first introduced by Yule (1926), and re-examined by Granger and Newbold (1974) using the Monte
Carlo simulations (Mills, 2019: p.235, 238).
53
where uYt and uXt are independent N(0,1) random errors.
Then the following OLS regression:
Yt =  + βXt + t (7.3)
Since Yt and Xt are independent of each other, any regression between them should give
insignificant results. However, after refreshing various Yts on Xts, they surprisingly
found that the null hypothesis of β2 = 0 was not rejected for approximately 75% of the
cases. In these simulations, characterized by a high R2 and a low value of Durbin-
Watson (dw) statistic, the usual t- and F-tests on the regression parameters may be
misleading (Verbeek, 2017: p.352).
To see the spurious regression problem, we can type the following commands in Stata:
clear
set obs 500
gen time = _n
set seed 12345
drawnorm e1 e2, n(500) means(0 0) sds(1 1)
tsset time
gen Y = 0
gen X = 0
replace Y = L.Y + e1 if _n > 1
replace X = L.X + e2 if _n > 1
label variable Y "Y is a random walk"
label variable X "X is a random walk"
twoway scatter Y X || lfit Y X, ytitle("Y is a random walk") xtitle("X is a random
walk") legend(off)
reg Y X
An example of a plot of Y against X obtained in this way is shown in Figure 7.1. The
estimated equation between these simulated series is presented in Table 7.1.
40
30
Y is a random walk

20
10
0

-20 -10 0 10 20
X is a random walk

Figure 7.1: Scatter plot of a spurious regression.

54
Table 7.1: A spurious regression with two random walks.

A similar example can be found in Hill et al. (2018: p.575), where t-statistic is huge
(40.84). These results are completely meaningless or spurious. The apparent
significance of the relationship is false. Hill et al. (2018: p.575) said that such spurious
results come from the fact that we have related one series with a stochastic trend to
another series with another stochastic trend.
Granger and Newbold (1974) proposed the following “rule of thumb” for detecting
spurious regressions: If R2 > DW statistic or if R2  1 then the estimated regression
‘must’ be spurious (Gujarati, 2011: p.226; Mills, 2019: p.237).
To understand the problem of spurious regression better, it might be useful to use an
example with real economic data. This example was conducted by Asteriou and Hall
(2016: p.353). Consider a regression of the logarithm of real GDP (Yt) to the logarithm
of real money supply (Mt) and a constant. The results obtained from such a regression
are the following:
Yt = 0.042 + 0.453Mt; R2 = 0.945; dw = 0.221 (7.4)
(t-stat.) (4.74) (8.57)
Here we see very high t-ratios, with coefficients that have the correct signs and more or
less plausible magnitudes. The coefficient of determination is very high (R2 = 0.945),
but there is a high degree of serial correlation (dw = 0.221). This shows evidence of the
possible existence of spurious regression. In fact, this regression is totally meaningless
because the money supply data are for the UK economy, while the GDP data are for the
US economy. Therefore, although there should not be any significant relationship, the
regression seems to fit the data well, and this happens because the variables used in the
example are trended (i.e., nonstationary). So, Asteriou and Hall (2016: p.353)
recommends that econometricians should be very careful when working with trended
variables. Similar regression results can be found if you try to regress real GDP of one
country on real money supply (or any macroeconomic variable) of another country. You
can see similar examples in Gujarati (2011, pp.224-226).

55
7.2 Explaining the Spurious Regression Problem
According to Asteriou and Hall (2016: p.353-4), if two variables, Xt and Yt, are both
stationary, then in general any linear combination of them will certainly be stationary.
One important linear combination of them is of course the error term, and so if both
variables are stationary, the error term will also be stationary and have a well-behaved
probability distribution (i.e., a white noise). However, when the variables are
nonstationary, then we cannot guarantee that the errors will be stationary and as a
general rule (although not always) the error itself is nonstationary. If this happens, we
violate the classical linear regression model (CLRM) assumptions of OLS regression
[e.g., zero mean value of error term, homoscedasticity, no serial correlation, that is, t
in equation (7.3) ~ idd(0,2)]48. If the errors are nonstationary, we could expect them to
wander around and eventually get larger. But OLS regression, it selects the parameters
so as to make the sum of the squared errors as small as possible, will select any
parameter which gives the smallest error and so almost any parameter value can result.
The simplest way to examine the behavior of t is to rewrite (7.3) as:
t = Yt –  – βXt (7.5)
or, excluding the constant  (which only affects t sequence by rescaling it):
t = Yt – βXt (7.6)
If Yt and Xt are generated by equations (7.1) and (7.2), then if we impose the initial
conditions Y0 = X0 = 0 we get that:
t = Y0 + uY1 + uY2 + … + uYt – β(X0 + uX1 + uX2 + … + uXt)
or
εt = ∑Ti=1 uYi – β ∑Ti=1 uXi (7.7)
From equation (7.7), we realize that the variance of the error term [var(t)] will tend to
become infinitely large as t increases (why?)49. Hence, the assumptions of the CLRM
are violated, and therefore, any t test, F test or R2 are unreliable.
In terms of equation (7.3), there are four different cases to discuss (Asteriou and Hall,
2016: p.355; Hill et al., 2017: p.587):
Case 1: If both Yt and Xt are stationary [i.e., I(0)]50, we can estimate a regression
relationship between the levels of those variables without fear of

48
For more details, see Gujarati and Porter (2009: p.189).
49
Because both Yt and Xt are random walks without drift. If Yt and Xt are either random walks with drift or with
drift plus a time trend, then both mean and variance of εt will tend to become infitely large as t increases.
50
The series Yt is said to be integrated of order one, denoted I(1), if it becomes stationary after taking a first
difference. A nonstationary series is integrated of order d, denoted I(d), if it becomes stationary after being first
differenced d times (Greene, 2018: p.1023). A stationary series such as AR(1) and MA(1) is denoted as I(0). In
general, the order of integration of a series is the minimum number of times it must be differenced to make it
stationary (Hill et al., 2017: p.581). According Greene (2018: p.1024), An I(1) will typically be constantly growing
or wandering about with no tendency to revert to a fixed mean. Most macroeconomic flows and stocks that relate
to population size, such as output or employment, are I(1). An I(2) series is growing at an ever-increasing rate,
such as the price level. Series that are I(3) or greater are extremely unusual, but the do exist.

56
 encountering a spurious regression. In other words, the CLRM is appropriate,
and OLS estimates are BLUE. This case is popularly discussed in
econometrics courses for undergraduates (i.e., serial correlation topic).
However, stationarity at levels rarely happens in macroeconomic time series.
Therefore, this case is beyond the scope of this notes.
Case 2: Yt and Xt are integrated of the same order [often I(1)] and the t sequence
contains a stochastic trend (i.e., a random walk). In this case, we have
spurious regressions and it is often recommended to either re-estimate the
regression equation in the first differences. This will be discussed more
clearly in Section 10.6.
Case 3: Yt and Xt are integrated of the same order [often I(1)] and the t sequence is
stationary. In this special case, Yt and Xt are said to be cointegrated. The
concept of cointegration and error correction model will be discussed in
Sections 9.2, 10.4, and 10.5.
Case 4: Yt and Xt are integrated of different orders [often I(0) and I(1), or sometimes
I(0) to I(d)]. In this case, the regression equations may be meaningless.
However, we can specify the model in different approaches or select other
proxies to direct the analysis toward either case 2 or case 3. This will be
discussed in Sections 10.6.

57
 8.
TESTING FOR NONSTATIONARITY

8.1 Graphical Analysis

According to Gujarati and Porter (2009: p.749), before implementing formal tests, it is
always advisable to plot the time series under study. Such a plot (line graph of the level
and the first difference) gives an initial clue about the likely nature of the time series.
Such an intuitive feel is the starting point for formal tests of stationarity (i.e., choose the
appropriate test equation). Note that we should declare the time series data by using
tsset command before plotting the data, using tsline command. For examples:

A white note series

By using the following Stata’s commands:

clear
set obs 1000
set seed 12345
drawnorm e, n(1000) means(0) sds(1)
label variable e "A white noise series"
gen date = date("9/8/17","MDY", 2017)+_n-1
format date %td
tsset date
tsline e
4 2
A white noise series
-2 0-4

01jul2017 01jul2018 01jul2019 01jul2020

date

Figure 8.1: A white noise series.

58
A stationary series

A stationary around a zero-mean series

By using the following Stata’s commands:

clear
set obs 500
set seed 12345
drawnorm e, n(500) means(0) sds(1)
gen month = m(1976m9)+_n-1
format month %tm
tsset month
gen stationary = 0
replace stationary = 0.7*[Link] + e if _n > 1
label variable stationary "stationary around a zero mean"
egen zero = mean(stationary)
label variable zero "mean"
tsline stationary zero
6
4
2
0
-2
-4

1980m1 1990m1 2000m1 2010m1 2020m1

month

stationary around a zero mean mean

Figure 8.2: A “stationary around a zero mean” series (Yt = 0.7Yt-1 + ut).

59
A stationary around a nonzero mean series

By using the following Stata’s commands:

clear
set obs 500
set seed 12345
drawnorm e, n(500) means(0) sds(1)
gen month = m(1976m9)+_n-1
format month %tm
tsset month
gen stationary = 0
replace stationary = 1 + 0.7*[Link] + e if _n > 1
label variable stationary "stationary around a nonzero mean"
egen nonzero = mean(stationary)
label variable nonzero "mean"
tsline stationary nonzero
6
4
2
0
-2
-4

1980m1 1990m1 2000m1 2010m1 2020m1

month

stationary around a nonzero mean mean

Figure 8.3: A “stationary around a nonzero mean” series (Yt = 1 + 0.68Yt-1 + ut).

60
A nonstationary series

A trend stationary series

By using the following Stata’s commands:

clear
set obs 500
set seed 12345
drawnorm e, n(500) means(0) sds(1)
gen time = _n-1
tsset time
gen Y = 0
replace Y = 1.2 + 0.01*time +0.7*L.Y + e if _n > 1
label variable Y "a trend stationary series"
reg Y time
predict trend
label variable trend "linear trend"
tsline Y trend
25
20
15
10
5
0

0 100 200 300 400 500

time

a trend stationary series linear trend

Figure 8.4: A trend stationary series (Yt = 1.2 + 0.01t + 0.7Yt-1 + ut).

61
A difference stationary series

By using the following Stata’s commands:

clear
set obs 500
set seed 12345
drawnorm e, n(500) means(0) sds(1)
gen time = _n
tsset time
gen Y = 0
replace Y = L.Y + e if _n > 1
label variable Y "a difference stationary series"
tsline Y
40
a difference stationary series
10 200 30

0 100 200 300 400 500

time

Figure 8.5: A nonstationary series (Yt = Yt-1 + ut).

By using the following Stata’s commands:
clear
set obs 500
set seed 12345
drawnorm e, n(500) means(0) sds(1)
gen time = _n
tsset time
gen Y = 0
replace Y = 0.1 + L.Y + e if _n > 1
label variable Y "a difference stationary series"
tsline Y

62
We have a time series plot of a difference stationary series:
80
a difference stationary series
20 400 60

0 100 200 300 400 500

time

Figure 8.6: A nonstationary series (Yt = 0.1 + Yt-1 + ut).

A seasonal series
3000000
2000000
sales
1000000
0

2005q1 2010q1 2015q1 2020q1

quarter

Figure 8.7: A seasonal series.

63
8.2 Tests for White Noise
In order to test whether a series (or a residual of a regression) is a white noise or not,
we may have several testing methods such as Bartlett’s periodogram-based test,
Portmanteau (Q) test.

Bartlett’s periodogram-based test for white noise

Bartlett’s test is a test of the null hypothesis that the data come from a white noise
process of uncorrelated random variables having a constant mean and a constant
variance. In Stata, we use the command: wntestb (i.e., wn = white noise, test = test, and
b = Bartlett).
clear
set seed 12393
set obs 100
gen time = _n
tsset time
gen x = rnormal()
gen y = 0
replace y = -0.6*L.y+rnormal() if _n>1
wntestb x
wntestb y
where x and y are the random and AR(1) series, respectively.

Cumulative Periodogram White-Noise Test

1.00 0.80
Cumulative periodogram for x
0.20 0.40 0.60
0.00

0.00 0.10 0.20 0.30 0.40 0.50

Frequency
Bartlett's (B) statistic = 0.80 Prob > B = 0.5380

Figure 8.8: A white noise series Xt.

64
We can see that the values in Figure 8.8 never appear outside the confidence bands. The
test statistic has a p-value (0.538), so we conclude that the process is not different from
white noise. If we had wanted only the statistic without the plot, we would have used
the table option.

Cumulative Periodogram White-Noise Test

1.00 0.80
Cumulative periodogram for y
0.20 0.40 0.60
0.00

0.00 0.10 0.20 0.30 0.40 0.50

Frequency
Bartlett's (B) statistic = 3.06 Prob > B = 0.0000

Figure 8.9: A random walk series, Zt.

In contrast, the values in Figure 8.9 appear outside of the bands. We also have
confirmation from the test statistic, at a p-value of 1 % level, that the Yt process is
statistically different from white noise.

Portmanteau (Q) test for white noise

Box and Pierce (1970) developed a portmanteau test of white noise that was refined by
Ljung and Box (1978). The statistic formulas are presented in Section 8.4. Note that the
Q statistic is also available with the corrgram command.
In Stata, we use the command: wntestq (i.e., wn = white noise, test = test, and q = Q).
The syntax is:
wntestq varname, lags(#)
where lags(#) specifies the number of autocorrelation to calculate. The default is to use
min ([n/2] – 2, 40), where [n/2] is the greatest integer less than or equal to n/2.
clear
set seed 12393
set obs 100
65
 gen time = _n
tsset time
gen x = rnormal()
gen y = 0
replace y = -0.6*L.y+rnormal() if _n>1
wntestb x
wntestq x
wntestq x, lags(10)

Both Bartlett and portmanteau tests result in nonsignificant test statistics: a p-value of
0.5380 for wntestb and a p-value of 0.8443 for wntestq. Therefore, Xt is a white noise.
wntestb y
wntestq y
wntestq y, lags(10)

Both Bartlett and portmanteau tests result in significant test statistics: a p-value of
0.0000 for both tests. This indicates that we reject the null hypothesis of white noise at
even 0.01% level.
66
8.3 Autocorrelation Function and Correlogram
Autocorrelation is the correlation between a variable lagged one or more periods and
itself. The correlogram is a graph of the autocorrelations for various lags of a time series
data51. According to Hanke (2014: p.22), the autocorrelation coefficients (see Section
5.2) of a series can be used to answer the following questions52:
(1) Are the data random? (This is usually used for the diagnostic checking of the
residuals such as from ARIMA model or OLS regressions with stationary time
series data). However, Bartlett and portmanteau (Q) tests for white noise are
mostly used.
(2) Are the data stationary?
(3) Do the data have a trend (nonstationary)?
(4) Are the data seasonal?
Besides, the correlogram is very useful when selecting the appropriate lags [i.e., p and
q] in the ARIMA models53.
(1) If a series is random, the autocorrelations (i.e., ACF) between Yt and Yt-k for
any lag k are close to zero (i.e., individual autocorrelation coefficients are
statistically insignificant or inside the shaded area). The successive values of a
time series are not related to each other (Figure 8.10). In other words, Yt and
Yt-k are completely independent for all values of k (k = 1, …., p).
(2) If a series is stationary, the autocorrelation coefficients for, say lag 1, lag 2, or
lag 3, are significantly different from zero and then suddenly die out as the
number of lags increases (Figures 8.11 and 8.12). In other words, Yt and Yt-1,
Yt and Yt-2, Yt and Yt-3 are weakly correlated [i.e., weakly dependent]; but Yt
and Yt-k [as k increases] are completely independent.
(3) If a series has a (stochastic) trend, successive observations are highly
correlated, and the autocorrelation coefficients are typically significantly
different from zero for the first several time lags and then gradually drop
toward zero as the number of lags increases [i.e., not weakly dependent]. The
autocorrelation coefficient for lag 1 is often very large (close to 1). The
autocorrelation coefficient for lag 2 will also be large, and so on. However, it
will not be as large as for lag 1 (Figure 8.13).
(4) If a series has a seasonal pattern, a significant autocorrelation coefficient will
occur at the seasonal time lag or multiples of seasonal lag (Figure 8.14). This
is beyond the scope of this series of lectures.
[

51
Correlogram also includes partial autocorrelation function.
52
For details, please see Hanke and Wichern (2014: p. 22-30).
53
As discussed in Section 4 about AR(p) and MA(q) that p is selected by using PACF graph, and q by using ACF
graph. ARIMA(p,d,q) is just a combination of the two processes for the series that is integrated of order d (see
Appendix B).
67
Figure 8.10: Correlogram of a purely random series54.

Figure 8.11: Correlogram of a stationary series: AR(1).

54
The Stata’s commands:
clear
set obs 10000000
set seed 12345
drawnorm whitenoise, n(10000000) means(0) sds(1)
gen time = _n
tsset time
label variable whitenoise "A white noise series"
corrgram whitenoise, lag(12)
68
Figure 8.12: Correlogram of a stationary series: MA(1).

Figure 8.13: Correlogram of a nonstationary series.

Figure 8.14: Correlogram of a seasonal series.

69
The correlogram becomes very useful for time series forecasting and other practical
(business) implications. For examples, after estimating ARIMA models, we often use
correlogram in order to check whether the error term is white noise. If it is not random
[e.g., the fourth AC coefficient (of the residuals) is statistically different from zero], we
have to re-estimate the model that includes the MA(4) term to account for quarterly
effect if the series was observed quarterly. If we conduct academic studies, however, it
is necessary to provide more formal statistics such as Z statistic, Box-Pierce Q statistic,
Ljung-Box (Q) portmanteau test, and especially unit root tests.

8.4 Tests for Nonstationarity

Barlett test for significance of an individual autocorrelation coefficient

According to Hill et al. (2018: p.425), it is often useful to test whether a sample
autocorrelation ρ̂𝑘 (see Section 5.2) is significantly different from zero. That is, a test of
H0: k = 0 against the alternative H1: k  0. Tests of this nature are useful for
constructing models (such as ARIMA models) and for checking whether the errors in
an equation might be serially correlated (or randomly distributed). Barlett (1946) stated
that if a time series is purely random, that is, it exhibits white noise, the sample
autocorrelation coefficients ρ̂𝑘 are approximately:

ρ̂𝑘 ~ N(0, 1/T) (8.1)

that is, in large samples the sample autocorrelation coefficients are normally distributed
with zero mean and variance equal to one over the sample size (Gujarati, 2009: p.753).
Therefore, the standard error of sample autocorrelation coefficients is simply √1⁄T.
According to Hill et al. (2018: p.426), the test statistic for the null hypothesis is:
̂k −0
ρ
Z= = ρ̂𝑘 √T ~ N(0, 1) (8.2)
√1/T

Following the properties of the standard normal distribution, the 95 percent confidence
interval for any (population) k is:

–1.96  ρ̂𝑘 √T  1.96 (8.3)

In other words,

Prob(ρ̂𝑘 – 1.96/√T  k  ρ̂𝑘 + 1.96/√T) = 0.95 (8.4)

According to equation (8.3), at 5% significance level, we reject H0: k = 0 when ρ̂𝑘 √T

 –1.96 or ρ̂𝑘 √T  1.96 (i.e., outside the 95% confidence interval). According to
equation (8.4), if this interval includes the value of zero, we do not reject H0: k = 0, but
if this interval does not include 0, we reject H0: k = 0. For example, suppose that total
observation is 244 and we test whether the coefficient of 0.78 at lag 10 is statistically
significant? Since we have 244 observations, the SQRT(244) = 15.62 and the product

70
ρ̂𝑘 √T = 0.78*15.62 = 12.18 > 1.96. Applying the equation (8.4), the 95% confidence
interval for true 10 is (0.62, 0.94). Obviously, this interval does not include the value
of zero. Therefore, we can reject H0: k = 0.
In Stata, we use the command gen z = sqrt(e(N))rk or gen z = sqrt(r(N))rk to generate
the confidence interval. For example,
import excel using D:\Table13_6.xls, firstrow clear
tsset time
ac close, lags(12) gen(ac_close)
qui sum close
qui return list
gen z = sqrt(r(N))*ac_close
gen lower = ac_close-1.96/sqrt(r(N))
gen upper = ac_close+1.96/sqrt(r(N))
list ac_close z lower upper in 1/4

Tests for significance of a group of autocorrelation coefficients

Instead of testing the statistical significance of any individual autocorrelation

coefficient, we can test the joint hypothesis that all the k up to certain lags are
simultaneously equal to zero (Gujarati, 2009: p.753). This can be done by using the Q
statistic developed by Box and Pierce (1970), which is defined as:

Q = T ∑m ̂ 2k
k=1 ρ (8.5)
where T = sample size and m = lag length. The Q statistic is often used as a test of
whether a time series is white noise (Gujarati, 2009: p.754). In order words, this is often
used for diagnostic checking of time series model specification such as ARIMA models,
that is, it tests whether the residual or forecast error is purely random. If Q statistic is
applied for testing stationarity, we have to consider the first lag lengths. According to
Gujarati (2009: p.754), in large samples, the Q statistic is approximately distributed as
the chi-square distribution with m df (i.e., the degree of freedom). If the computed Q
exceeds the critical Q value from the chi-square distribution at the chosen level of
significance, one can reject the null hypothesis that all the (true) k are zero; at least
some of them must be nonzero.

71
A variant of the Box-Pierce Q statistic is the Ljung-Box (LB) statistic (Ljung and Box,
1978), which is defined as:
̂2k
ρ
LB = T(T + 2)∑m
k=1 (
2
) ~ 𝜒m (8.6)
T−k

Although in large samples both Q and LB statistics follow the chi-square distribution
with m df, the LB statistic has been found to have better (more powerful, in the statistical
sense) small-sample properties than the Q statistic (Gujarati, 2009: p.754).
The Stata’s command for Q statistic is the corrgram, while the Stata’s command for LB
statistic is the lmalb55. For example,
import excel using D:\Table13_6.xls, firstrow clear
tsset time
corrgram [Link], lag(4)
lmalb [Link], lag(4)

55
We also use the Ljung-Box (Q) portmanteu test (see Appendix B).

72
Simple Dickey-Fuller test for unit roots

A test for non-stationarity that has become widely popular over the past several decades
is the unit root test (Gujarati, 2009: p.754). Dickey and Fuller (1979, 1981) proposed a
procedure to formally test for non-stationarity (hereafter refer to Dickey-Fuller or DF
tests). The key insight of these tests is that testing for non-stationarity is equivalent to
testing for the existence of a unit root.
There are three variations of the Dickey-Fuller tests, each one designed for a different
alternative hypothesis (Hill et al., 2018: p.577):
(1) The alternative hypothesis is that Yt is stationary around a zero mean as
illustrated in Figure 8.2. Both intercept and trend are excluded from the test
equation in this case.
(2) The alternative hypothesis is that Yt is stationary around a nonzero mean as
illustrated in Figure 8.3. In this case, the test equation includes an intercept but
no trend term.
(3) The alternative hypothesis is that Yt is stationary around a linear trend as
illustrated in Figure 8.4. The test equation includes both intercept and trend term
in this case. In other words, the alternative hypothesis is that Yt is trend
stationary.

Dickey-Fuller test without intercept and trend

Suppose we wish to test whether a time series Yt is better represented by a stationary

AR(1) process like that in Figure 8.2 or a random walk process like that in Figure 8.5.
The nonstationary random walk (i.e., a random walk without drift) is set up as the null
hypothesis:
H0: Yt = Yt-1 + ut (8.7)
and the stationary AR(1) process becomes the alternative hypothesis:
H1: Yt = 1Yt-1 + ut, where |1| < 1 (8.8)
Throughout, we assume that the ut are independent random errors with mean zero and
variance 2, and they are uncorrelated with the past values Yt-1, Yt-2, , …, Yt-p. That is,
ut ~ N(0, 2) and Cov(ut, ut) = 0  t  s (Greene, 2018: p.1029). Under H1, the series
fluctuates around zero mean. Under H0, it wanders upward and downward but does not
exhibit a clear trend in either direction and does not tend to return to a constant mean.
However, an obvious way to specify the hypotheses in terms of the parameters such as
H0: 1 = 1 against H1: |1| < 1 or simply H1: 1 < 1. Therefore, equation (8.8) is rewritten
as follows:
Yt = 1Yt-1 + ut, where |1|  1 (8.9)
However, we cannot estimate equation (8.9) by OLS and test the null hypothesis that 
= 1 by the usual t test because such a test is severely biased due to spurious regression
in the case of a unit root (Gujarati, 2009: p.754). Therefore, we have to manipulate
equation (8.9) as follows: Subtract Yt-1 from both sides of equation (8.9) to obtain:
73
 Yt – Yt-1 = 1Yt-1 – Yt-1 + ut
∆Yt = (1 – 1)Yt-1 + ut
∆Yt = Yt-1 + ut (8.10)
where  = (1 – 1) and ∆Yt = Yt – Yt-1. Therefore, instead of estimating equation (8.9),
we estimate equation (8.10) and test the null hypothesis that H0:  = 0 against the
alternative hypothesis that H1:  < 0. Rejection of the null hypothesis that  = 0 implies
the series is stationary around a zero mean. A failure to reject H0 suggests the series
could be nonstationary, and we must be careful not to proceed to estimate a spurious
regression (Hill et al., 2018: p.578).
Now let us turn to the estimation of equation (8.10). This is rather simple; all we have
to do is to take the first difference of Yt, and regress it on Yt-1 and see if the estimated
slope coefficient in this regression (δ̂) is zero or not. If it is zero, we conclude that Yt is
nonstationary. But if it is negative, we conclude that Yt is stationary56.

Dickey-Fuller test with intercept and no trend

Suppose we wish to test whether a time series Yt is better represented by a stationary

AR(1) process like that in Figure 8.3 or a nonstationary random walk like that in Figure
8.6. The nonstationary random walk (i.e., a random walk with drift) is set up as the null
hypothesis:
H0: Yt =  + Yt-1 + ut (8.11)
and the stationary AR(1) process becomes the alternative hypothesis:
H1: Yt =  + 1Yt-1 + ut, where || < 1 (8.12)
If  = 0, equations (8.11) and (8.12) will become equations (8.7) and (8.8). Under H0
and H1 in terms of the parameters, equation (8.12) is rewritten as:
Yt =  + 1Yt-1 + ut, where |1|  1 (8.13)
or
∆Yt =  + Yt-1 + ut (8.14)
where  = (1 – 1) and ∆Yt = Yt – Yt-1. Therefore, instead of estimating equation (8.13),
we estimate equation (8.14) and test the null hypothesis that H0:  = 0 against the
alternative hypothesis that H1:  < 0. Rejection of the null hypothesis that  = 0 (and 
= 0) implies the series is stationary around a nonzero mean. A failure to reject H0
suggests the series could be nonstationary.
Now let us turn to the estimation of equation (8.14). This is rather simple; all we have
to do is to take the first difference of Yt, and regress it on Yt-1 and an intercept, and see
if the estimated slope coefficients in this regression (δ̂ and α ̂) are zero or not. If δ̂ is
zero, we conclude that Yt is nonstationary. But if it is negative and α ̂ is statistically
different from zero, we conclude that Yt is stationary around a nonzero mean.

56
Since  = (1 - 1), for stationarity  must be less than one. For this to happen  must be negative.
74
Dickey-Fuller test with intercept and trend

Above tests are used if a time series Yt has a trend upward or downward. In case, as
illustrated in Figure 8.4, Yt might be stationary around a linear trend and described by
the process:
Yt =  + t + 1Yt-1 + ut, where || < 1 (8.15)
where t is the time or trend variable. A time series that can be described by equation
(8.15) is called trend stationary if  is different from zero. We already know that a
random walk with drift, as illustrated in Figure 8.7, is described as:
Yt =  + Yt-1 + ut (8.16)
In this case, Yt is nonstationary. The Dickey-Fuller test with intercept and trend is
designed to discriminate between these two models. Equation (8.16) becomes the null
hypothesis (H0), and equation (8.15) is the alternative hypothesis (H1). If the null
hypothesis is rejected, we conclude Yt is trend stationary. Failure to reject H0 suggests
Yt is nonstationary, or at least there is insufficient evidence to prove otherwise (Hill et
al., 2018: p.579).
If  = 0, equations (8.15) and (8.16) will become equations (8.11) and (8.12). Under H0
and H1 in terms of the parameters, equation (8.14) is rewritten as:
Yt =  + t + 1Yt-1 + ut, where ||  1 (8.17)
or
∆Yt =  + t + Yt-1 + ut (8.18)
where  = (1 – 1) and ∆Yt = Yt – Yt-1. Therefore, instead of estimating equation (8.17),
we estimate equation (8.18) and test the null hypothesis that H0:  = 0 against the
alternative hypothesis that H1:  < 0. Rejection of the null hypothesis that  = 0 (and 
= 0) implies the series is stationary around a linear trend. A failure to reject H0 suggests
the series could be nonstationary.
Now let us turn to the estimation of equation (8.18). This is rather simple; all we have
to do is to take the first difference of Yt, and regress it on Yt-1, an intercept and a time
variable, and see if the estimated slope coefficients in this regression (δ̂ and γ̂) are zero
or not. If δ̂ is zero, we conclude that Yt is nonstationary. But if it is negative and γ̂ is
statistically different from zero, we conclude that Yt is stationary around a linear trend.
In short, if the null hypothesis is rejected, it means either: (1) Yt is stationary with zero
mean, in the case of equation (8.7), or (2) Yt is stationary with nonzero mean, in the
case of equation (8.11). In the case of equation (8.18), we can test for  < 0 (i.e., no
stochastic trend) and   0 (i.e., the existence of a deterministic trend) simultaneously,
using the F test, but using the critical values tabulated by Dickey and Fuller (Gujarati,
2009: p.756).
According to Gujarati (2009: p.755), the only question is which test we use to find out
if the estimated coefficient of Yt-1, that is δ̂, in equations (8.10, 8.14 and 8.18) is zero
or not. Under the null hypothesis that H0:  = 0 (i.e., 1 = 1), the t value of the estimated
coefficient of Yt-1 no longer follows the t distribution even in large samples. According

75
to Hill et al. (2018 : p.578), when H0 is true, Yt is not stationary and has a variance that
increases as the sample size increases. This increasing variance alters the distribution
of the usual t-statistic when H0 is true. Therefore, we cannot use the usual t test.
What is the alternative? Dickey and Fuller (1979) have shown that under the null
hypothesis that  = 0, the estimated t value of the coefficient of Yt-1 in equations (8.10,
8.14, and 8.18) follows the  (tau) statistic. These authors have computed the critical
values of the tau statistic on the basis of Monte Carlo simulations. MacKinnon (1991,
1996) prepared more extensive tables, which are now incorporated in several statistical
software packages. The tau statistic is widely known as the Dickey-Fuller (DF) test.
Keep in mind that the Dickey-Fuller test is one-sided because the alternative hypothesis
is that  < 0 (or 1 < 1). If the hypothesis that  = 0 is rejected (i.e., the time series is
stationary), we can use the usual t test (Gujarati, 2009 : p.755). Some critical values for
each of the preceding three specifications of the DF test are presented in Table 8.1.
Table 8.1: Critical values for DF tests with large sample size57.
Model 1% 5% 10%
∆Yt = Yt-1 + ut –2.58 –1.95 –1.62
∆Yt =  + Yt-1 + ut –3.42 –2.86 –2.57
∆Yt =  + t + Yt-1 + ut –3.96 –3.41 –3.13
Standard normal critical values –2.33 –1.65 –1.28
Source: Asteriou & Hall (2016: p.356), Hill et al. (2018: p.578)58.
It is extremely important to note that the critical values of the tau statistic to test the
hypothesis that  = 0 is different for each of the preceding three specifications of the
DF test. There is no way of knowing which specification is correct to begin with. Some
trial and error is inevitable, data mining notwithstanding (Gujarati, 2009 : p.756). The
actual estimation procedure is as follows: (i) Estimate equation (8.10), or equation
(8.14), or equation (8.18) by OLS; (ii) divide the estimated coefficient of Yt-1 in each
case by its standard error to compute the tau statistic; and (iii) refer to the DF table (or
any statistical package).
If the computed absolute value of the tau statistic (| |) exceeds the absolute DF or
MacKinnon critical tau values, we reject the hypothesis that  = 0, in which case the
time series is stationary. On the other hand, if the computed | | does not exceed the
absolute critical tau value, we do not reject the null hypothesis, in which case the time
series is nonstationary. In most applications the tau value will be negative. Therefore,
alternatively we can say that if the computed (negative) tau value is smaller than (i.e.,
more negative than) the critical tau value, we reject the null hypothesis (i.e., the time
series is stationary), otherwise, we do not reject it (i.e., the time series is nonstationary).
An easy way is to compare the ‘MacKinnon approximate’ p-value and the significance
level () often at 1%, 5%, and 10%. If the p-value is smaller than a chosen level of
significance, we reject the null hypothesis of a unit root. Note that the MacKinnon

57
For more details, see Greene (2018: p.1030).
58
For more details, see Fuller (1996: p.642).
76
approximate p-value and the test statistic are not always consistent to each other (see
StataCorp, 2019a: dfgls).

Augmented Dickey-Fuller test for unit roots

In conducting the DF test as in equations (8.10), (8.14), and (8.18), it is assumed that
the error term ut is uncorrelated. But in case the ut are correlated (i.e., serial correlation),
Dickey and Fuller have developed another test, known as the augmented Dickey-Fuller
(ADF) test. This test is conducted by “augmenting” the preceding three equations by
adding the lagged values of the dependent variable Yt (Gujarati, 2009: p.757). In other
word, the ADF test is used in order to deal with serial correlation in the test equations59.
The ADF test has three possible specifications:
p
∆Yt = Yt-1 + ∑i=1 β𝑖 ∆Yt−i + ut (8.19)
p
∆Yt =  + Yt-1 + ∑i=1 β𝑖 ∆Yt−i + ut (8.20)
p
∆Yt =  + t + Yt-1 + ∑i=1 β𝑖 ∆Yt−i + ut (8.21)
where ut is a pure white noise error term and where Yt = Yt – Yt-1, Yt-1 = Yt-1 – Yt-2,
Yt-2 = Yt-2 – Yt-3, etc.
Similar to the DF cases, the ADF tests concern the null hypothesis that H0:  = 0 after
the test equation is correctly specified. Keep in mind that if the computed absolute value
of the tau statistic (| |) exceeds the absolute ADF or MacKinnon critical tau values, we
reject the hypothesis that  = 0, in which case the series is stationary. On the other hand,
if the computed | | does not exceed the absolute critical tau value, we do not reject the
null hypothesis, in which case the time series is nonstationary. In most applications the
tau value will be negative. Therefore, alternatively we can say that if the computed
(negative) tau value is smaller than (i.e., more negative than) the critical tau value, we
reject the null hypothesis (i.e., the time series is stationary), othereise, we do not reject
it (i.e., the time series is nonstationary).

Selection of optimal lag length

The inference of ADF test is very sensitive to the selection of lag length (Schwert, 2002;
Agiakloglou and Newbold, 1992; Harris, 1992). In carrying out ADF test, we have to
determine the optimal lag length p which depends on the data generating process of the
series. It is shown by Ng and Perron (1995) that if p is chosen on the basis of some
information criteria (IC) and the general special criteria (GSC) the asymptotic
distribution of tδ will follow the standard DF distribution. The method of selecting lag
length on the basis of information criteria considers a trade-off between the size
distortions because of the inclusion of too few lags and the power losses because of the
inclusion of too many lags (Das, 2019: p.320).
The optimal lag length (p*) is obtained as:

59
More exactly, the ADF is used if the test equation follows the AR models with higher order such as AR(2),
AR(3), …, AR(p) rather than the AR(1) model A(see Appendix D for more details).
77
 p* = arg min IC(p) (8.22)
pmin  p  pmax
CT
̂2p + p
IC(p) = ln σ (8.23)
T

̂2p is the OLS residuals from the pth order ADF regression, CT is a penalty function
Here σ
defined differently in different information criteria to be used. The most popularly used
information criteria are the Akaike information criterion (AIC) proposed by Akaike
(1974) and the Bayesian information criterion (BIC) proposed by Schwartz (1978). In
AIC, CT = 2, while for BIC, CT = lnT.
̂2p + 2p
AIC = Tln σ (8.24)
̂2p + pln(T)
BIC = Tln σ (8.25)
For finite p, the properties of AIC and BIC remain the same when Y t is stationary. In
presence of unit roots, however, BIC is consistent while AIC is not (Tsay, 1984).
Ng and Perron (1995) suggest the procedure for selecting lag length in such a way that
results in minimal power loss. In this procedure, we have to set an upper bound pmax for
p and estimate the ADF regression with p = pmax. If the coefficients for last lagged
differences are not statistically significant, we need to reduce the lag length one by one
and repeat the process. A useful rule of thumb for determining p max, suggested by
Schwert (2002), is:
T 0.25
pmax = [ 12 x ( ) ] (8.26)
100

Ng and Perron (2001) propose the following class of modified information criteria
(MIC):
CT+τT (𝑝)
̂2p + p
MIC(p) = ln σ (8.27)
T

where
T
2
̂2p ) ∑ Yt−1
τT (p) = 1/( σ (8.28)
pmax +1

The penalty function CT = 2 yields the modified AIC (MAIC), and CT = lnT yields the
modified BIC (MBIC). The value of τT(p) will decrease the increase of p.
The number of lagged difference terms to include is sometimes determined
empirically60, the idea being to include enough terms so that the error term in above
equations is serially uncorrelated (Gurajati, 2009: p.757). Stata has an option that
automatically selects the lag length based on Akaike, Schwarz, and other information
criterion. After each estimation, we may check whether the residuals of the ADF

60
We may also use the PAC graph of Yt to identify the possible lags of the test equation. For examples, if only
PAC1 and PAC2 are statistically different from zero, we may add only one lagged difference term; if PAC1, PAC2
and PAC3 are statistically different from zero, we may add two lagged difference terms; etc (see Appendix D).
An useful example of lag length selection is presented in Appendix B.
78
regression are serially correlated or not through either LM test61, Q test, or LB test. In
addition, we can also use the PAC graph to check whether any coefficient is still
different from zero.
In Stata, varsoc is often used to find out the optimal lag length. However, the pre-
estimation version of varsoc requires a maximum lag order. The command varsoc
provides the optimal lag length on the basis of final prediction error (FPE), Akaike
information criterion (AIC), Hannan and Quinn information criterion (HQIC), and
Schwartz Bayesian information criterion (SBIC). Although FPE is not an information
criterion, we want to minimise the prediction error to find out optimal leg length. It is
to be noted here that AIC is more accurate for monthly data, HQIC works better for
quarterly data on samples over 120, and SBIC works fine with any sample size for
quarterly data (Das, 2019: p.321). The command varsoc also provides the likelihood
ratio (LR) test statistics for the series of order less than or equal to the highest lag order
on the basis of log-likelihood (LL) function.
For univariate model, the LL at lag length p is:
̂2p
LL = σ (8.29)
For a given lag p, the LR test compare the AR(p) model with AR(p - 1) model:
LR(p) = 2{LL(p) – LL(p - 1)} (8.30)
The null hypothesis is that the coefficient on the pth lag of the endogenous variable is
zero.
Here FPE at lag length p is defined as:
T+p+1
̂2p (
FPE = σ ) (8.31)
T−p−1

The Stata commands are:

use "D:\Table13_1.dta", clear
tsset time
varsoc lnex, maxlag(6)
Table 8.2: Lag length selection in Stata.

61
Because Durbin-Watson d test cannot be used for higher order of autocorrelation (see basic econometrics
textbooks if you forget this point).
79
Figure 8.15: Procedure for testing for unit roots using ADF methodology.

Estimate the model

NO STOP: Conclude that

 = 0? there is no unit root

YES: Test for the presence of

the trend NO

is  = 0? given NO STOP: Conclude that

that  = 0?  = 0?
YES Yt has a unit root

YES

Estimate the model STOP: Conclude that

NO
there is no unit root

is  = 0?

YES: Test for the NO

presence of the constant

is  = 0? given NO STOP: Conclude that

that  = 0?  = 0? Yt has a unit root
YES

YES
STOP: Conclude that
there is no unit root
Estimate the model NO

is  = 0?
YES STOP: Conclude that
Yt has a unit root
Source: Asteriou and Hall (2016: p.358).

80
According to Asteriou and Hall (2016: p.357), unless the econometrician knows the
actual data-generating process, there is a question concerning whether it is most
appropriate to estimate models in (8.19), or (8.20), or (8.21). It is suggested that the test
procedure should start estimating the most general model given by equation (8.21) and
then answering a set of questions regarding the appropriateness of each model and
moving to the next model. This procedure is illustrated in Figure 8.15.
It needs to be stressed here that, although useful, this procedure is not designed to be
applied in a mechanical fashion. Plotting the data and observing the graph is sometimes
very useful because it can clearly indicate the presence or not of deterministic regressors
(StataCorp, 2019a: dfuller). [Note: we mean tsline of the first differenced series].
However, this procedure is the most sensible way to test for unit roots when the form
of the actual data-generating process is typically unknown. In addition, the ADF test
results are sensitive to the lag lengths selected (StataCorp, 2019a: dfgls). Therefore, in
practical applications, it is necessary to use other tests for comparison purpose62.

DF-GLS test of unit roots

Several tests have explicitly concentrated on improving power (i.e., it is not always to
tell if a unit root exists because DF and ADF tests have low power against near unit-
root alternatives, e.g., 1 = 0.95 for AR(1), 1 + 2 = 0.95 for AR(2), 1 + 2 + 3 = 0.95
for AR(3), etc) and reducing size distortion (i.e., because we cannot include an infinite
number of augmentation lags in ADF test equation). These tests are based on
generalized least squares (GLS) “detrending” prior to calculating a test statistic. The
Dickey-Fuller using GLS (DF-GLS), developed by Elliott, Rothenberg, and Stock
(1996), is based on this approach (Mills, 2019: p.84). Among statistical tests for a unit
root, DF-GLS is the most powerful and informative (Hamilton, 2012: p.376; StataCorp,
2019a: dfgls). A test with higher power than the ADF test is more likely to reject the
null hypothesis of a unit root against the stationary alternative when the alternative is
true; thus a more powerful test is better able to distinguish between a unit AR root and
a root that is large but less than 1 (Stock and Watson, 2015: p.697)63. Basically, DF-
GLS is an augmented Dickey-Fuller test, except that the series is transformed via a
generalized least squares regression before performing the test (see StataCorp, 2019a:
dfgls). This test is introduced for the case that, under the null hypothesis, Yt has a
random walk trend, possibly with drift, and under the alternative hypothesis, Yt is
stationary around a linear time trend (Stock and Watson, 2015: p.697).
According to Stock and Watson (2015: p.699), the DF-GLS test is computed in two
steps. In the first step, the intercept and trend are estimated by generalized least squares.
The GLS estimation is performed by computing three new variables, Vt, X1t, and X2t in
the following manner:

62
See Appendix B: Example 3.
63
For example, a time series with ∅ ̂ = 0.95 or 𝛿̂ = 0.05 is in fact a stationary AR(1), but it is more likely to be
nonstationary if we use the ADF test. Stock and Watson (2015: p.699) implemented an AR(1) simulation with 
= 0.95, T = 200, and without a time trend. They found that the probability that the ADF test correctly rejects the
null hypothesis at the 5% significance level is approximately 31% compared to 75% for the DF-GLS test.
81
 Y1 if t = 1
Vt = { (8.32)
Yt − 𝛼 ∗ Yt−1 if t = 2, … , T
1 if t = 1
X1t = { (8.33)
1 − 𝛼∗ if t = 2, … , T
1 if t = 1
X2t = { (8.34)
t − 𝛼 ∗ (t − 1) if t = 2, … , T
where 𝛼 ∗ is computed using the formula 𝛼 ∗ = 1 – 13.5/T. Then Vt is regressed against
X1t and X2t; that is, OLS is used to estimate the coefficients of the following population
regression equation:
Vt = 1X1t + 2X2t + et (8.35)
using the observations t = 1, 2, …, T, where et is the error term. Note that there is no
intercept in the regression in equation (8.35). The OLS estimators φ̂ 1 and φ
̂ 2 are then
d
used to compute a “detrended” version of Yt, namely, Yt :
Ytd = Yt – (φ
̂1 + φ
̂ 2 t) (8.36)
In the second step, the (augmented) Dickey-Fuller test is used to test for a unit
autoregressive root in Ytd , where the (augmented) Dickey-Fuller regression does not
include an intercept or a time trend, that is, we use either equation (8.10) or equation
(8.19). Note that the number of lags p is automatically determined by Stata by using
information criteria such as the MAIC, SBC.
Stock and Watson (2015: p.699) also notes that if the alternative hypothesis is that Yt is
stationary around a nonzero mean and without a time trend, the preceding steps are
modified. Specifically, 𝛼 ∗ is computed using the formula 𝛼 ∗ = 1 – 7/T, X2t is omitted
from the regression in equation (8.25), and the detrended series is computed as
Ytd = Yt – φ
̂1 (8.37)
Because the coefficients on the deterministic terms are estimated differently in the ADF
and DF-GLS tests, the tests have different critical values. Although the critical values
are routinely incorporated in statistical softwares such as Stata, here is an example:
Table 8.3: Critical values for DF-GLS test.
Model 1% 5% 10%
Intercept only (X1t only) –2.58 –1.95 –1.62
Intercept and time trend (X1t and X2t) –3.48 –2.89 –2.57
Source: Stock and Watson (2015: p.700).

Phillips-Perron test of unit roots

In practical studies, researchers mostly use both the ADF and the Phillips-Perron (PP)
tests. Because the distribution theory that supporting the Dickey-Fuller tests is based on
the assumption of random error terms [iid(0,2)], when using the ADF methodology we
have to make sure that the error terms are uncorrelated and they really have a constant
variance. Phillips and Perron (1988) have suggested an alternative to the augmented
82
Dickey-Fuller tests that allow for fairly mild assumptions concerning the distribution of
errors (Asteriou and Hall, 2016: p.357-9; Verbeek, 2017: p.306). The regression for the
PP test is similar to DF equation (8.17) or equation (8.18).
Yt =  + t + 1Yt-1 + ut, where |1|  1 (8.17)
or
∆Yt =  + t + Yt-1 + ut (8.18)
Instead of adding additional lags of dependent variable on the right-hand side in the
regressions such as the equation 8.21) to obtain an error term that has no autocorrelation
(i.e., white noise), the idea here is to estimate an “unaugmented” (or original) model –
say, equation (8.18) and then to adjust the DF test statistics to take into account the
(potential) autocorrelation pattern in the errors (Mills, 2019: p.83; Verbeek, 2017:
p.306). These adjustments, based on corrections similar to those applied to compute
Newey-West (HAC) standard errors64 to account for autocorrelation (Asteriou and Hall,
2016: p.359; Verbeek, 2017: p. 306; StataCorp, 2019a: pperron).
So, the PP statistics are just modifications of the ADF t statistics that take into account
the less restrictive nature of the error process. The expressions are extremely complex
to derive and are beyond the scope of my notes. Since most statistical packages have
routines available to calculate these statistics, it is good for researcher to test the order
of integration of a series performing the PP test as well. The asymptotic distribution of
the PP t statistic is the same as the ADF t statistic and therefore the MacKinnon
(1991,1996) critical values are still applicable. That means the PP tests also concern
whether  = 0. The PP test statistic is the t statistic for the lagged dependent variable. If
the PP statistical value is smaller [in absolute terms] than the critical value then we
reject the null hypothesis of a unit root and conclude that Yt is a stationary process.
Again, an easy way is to compare the MacKinnon approximate p-value and the
significance level () often at 1%, 5%, and 10%. If the MacKinnon approximate p-
value is smaller than a chosen level of significance, we reject the null hypothesis that
Yt represents a random walk or has a unit root. It is also important to know that the PP
tests can be performed with the inclusion of a constant and linear trend, or none of them
in the test regressions.

KPSS test of stationarity

To circumvent the problem that (augmented) Dickey-Fuller tests often have low power
(especially the series has a deterministic trend), Kwiatkowski, Phillips, Schmidt and
Shin (1992) propose an alternative test where stationarity is the null hypothesis and the
existence of a unit root is the alternative hypothesis. This test is usually referred to as
the KPSS test. According to Verbeek (2017: p.304), the basic idea is that a time series
is decomposed into the sum of a deterministic time trend, a random walk and a
stationary error term (typically not white noise). The null hypothesis (of trend
stationarity) specifies that the variance of the random walk component is zero. The test
is actually a Lagrange multiplier test, and computation of the test statistic is fairly

64
These are quite complicated and will not be discussed here. If you are interested, you can see them in
Wooldridge (2013: p.431-4) or Das (2019: p.327).
83
simple. According to Das (2019: p.331-2), Greene (2018: p.1038-9), and Verbeek
(2017: p.304), the procedure of this test is as follows:

(1) Run an auxiliary regression of Yt upon an intercept and a time trend t [i.e., Yt =
 + t + t, where t is a trend variable and t is an error term];
(2) Save the OLS residuals et and compute the partial sums St = ∑ts=1 es for all t;
(3) Calculate the test statistic
S2
KPSS = T −2 ∑T
s=1
t
(8.38)
̂2
σ
where σ ̂2 is an estimator for the long-run variance ∑∞ j=−∞ E{εt εt−j }. This
estimator is a weighted average of the sample autocovariance and several
alternative weighting schemes have been proposed. Most popular are the Barlett
weights used by KPSS and the quadratic spectral kernel. In practice, the KPSS
̂2 . (Verbeek,
test appears to be quite sensitive to the choices made to estimate σ
∑T 2
t=1 et
̂2 is given by σ
2017: p.304). According to Greene (2018: p.1039), σ ̂2 = +
T
j ∑T
s=j+1 es es−j
2 ∑Lj=1 (1 − ) rj where rj = and L is chosen by the analyst.
L+1 T
(4) Compare the test statistic with critical value at 5% significance level.
The critical values are routinely produced by statistical packages such as Stata. The null
hypothesis (of trend stationarity) is rejected if the KPSS test statistic is larger than the
selected critical value, often at 5% level of significance.

Zivot and Andrews test of unit roots (with endogenous break)

A problem common with the conventional unit root tests such as the ADF, DF-GLS, PP
and KPSS tests is that they do not account for possible structural breaks65 in estimation
of unit roots. As a result, they tend to bias the conclusion towards a unit root (i.e., in
favour of the null hypothesis of unit root) if there is a break in the series (Ling et al.,
2013; Arestis and Mariscal, 1999; Chaudhuri and Wu, 2003; Perron, 1989). The
structural break is assumed to be either exogenous66 (Perron, 1989) or endogenous67
(Zivot and Andrews, 1992) in the series68. However, the selection of an exogenous
structural break point will lead to over rejection of the unit root hypothesis, while the
inclusion of an unknown break in the series will reduce the bias in the unit root test
(Perron, 1997; Perron and Vogelsang, 1992; Zivot and Andrews, 1992). There are three
different models for testing unit roots with a single structural break:

65
Structural breaks are normally associated with anomalous events such as the oil crises in 1973-75, 1979-80, the
commodity crisis in 1985-86, the Asian financial crisis in 1997-98, the SARS epidemic in 2003, the energy shock
in 2005, the global financial crisis in 2008-09 (Ling et al., 2013).
66
This is also called ‘a known break’ (Mills, 2019: p.103) or ‘a predetermined break’ (Ling et al., 2013).
67
This is also called ‘an unknown break point’ (Mills, 2019: p.110) or ‘an estimated break’ (Ling et al., 2013).
68
For details, a good reference is Mills (2019: p.103-110).
84
 • Model A (model with intercept)69:
p
Yt =  + Yt-1 + t + DUt + ∑i=1 βi ∆Yt−i + ut (8.39)
• Model B (model with trend)70:
p
Yt =  + Yt-1 + t + DT*t + ∑i=1 βi ∆Yt−i + ut (8.40)
• Model C (model with both intercept and trend)71:
p
Yt =  + Yt-1 + t + DUt + DT*t + ∑i=1 βi ∆Yt−i + ut (8.41)
where the intercept dummy DUt = 1 if t > Tb and zero otherwise (where TB is the time
break), the slope dummy DT*t = t – Tb if t > TB and zero otherwise.
In equations from (8.39) to (8.41), the null hypothesis of a unit root (i.e., H0:  = 0)
against the alternative hypothesis that the series is stationary (i.e., Ha:  < 0). According
to Zivot and Andrews (1992), every point is considered as a potential break point (Tb)
and regression will be run for every possible break point sequentially. Among all the
possible break points, a break point (Tb) is selected when the absolute value of t-statistic
from the ADF test is minimized. In addition, the break point can also be searched over
the range of sample (0.15T – 0.85T).

Other unit roots with structural breaks

In this section, we introduce two additional tests of unit root that take structural breaks
into account. Perron and Vogelsang (1992) test takes one structural break in a series
into account. The Perron and Vogelsang (1992) unit root test has tow forms. The first
form of the test is the additive outlier model (AO) that captures a sudden change in the
series, if it exists. The second form is the innovative outlier (IO) model that captures
any gradual shift in the series. However, it has been argued that one endogenous break
is insufficient to identify the correct order of integration if the series contains more
breaks. The Perron and Vogelsang (1992) test fails to deal with more structural breaks.
Clemente, Montañés, and Reyes (1998) modifies the Perron and Vogelsang (1992) by
extending to two endogenous structural breaks in the series.

Perron and Vogelsang (1992) test

Additive outlier model:

For a fixed value of the breakpoint Tb, the test is a two-step procedure. First, we remove
the deterministic part of the series using the estimates of the regression:

Yt = μ + δDUt + ỹt , t = 1, 2, …, T (8.42)

Second, we perform the t statistic for  = 1 hypothesis in the following regression:

69
This is also called the ‘level shift’ model (Mills, 2019: p.103).
70
This is also called the ‘changing growth’ model (Mills, 2019: p.103).
71
This is also called the ‘combined’ model (Mills, 2019: p.104).
85
 𝑘 𝑘

ỹt = ∑ ωi D(TB)t−i + ỹt−1 + ∑ ci Δỹt−i + et

𝑖=0 𝑖=1

t = k + 2, …., T (8.43)

where DUt = 1 if t > Tb and 0 other wise, Tb is the break data, D(TB)t = 1 if t = Tb + 1
and 0 otherwise.
Innovational outlier model:
The IO model is specified as:
𝑘

Yt = μ + δDUt + θD(TB)t + Yt−1 + ∑ ci ΔYt−i + et (8.44)

𝑖=1

The IO model can be specified in other forms:

𝑘

Yt = μ + δDUt + βt + θD(TB)t + Yt−1 + ∑ ci ΔYt−i + et (8.45)

𝑖=1
𝑘

Yt = μ + δDUt + βt + γDTt + θD(TB)t + Yt−1 + ∑ ci ΔYt−i + et (8.46)

𝑖=1

where DTt = t – Tb if t > Tb and 0 otherwise.

Equations (8.45) and (8.46) are known as IO1 and IO2, respectively.

Clemente, Montañés, and Reyes (1998) test

This test objective is to extend the Perron and Vogelsang (1992) test to the case of two
changes in the mean of the variable being studied.
Additive outlier model:
The test procedure is as follows. First, we remove the deterministic part of the series by
estimating the following model:

Yt = μ + δ1 DU1t + δ2 DU2t + ỹt , t = 1, 2, …, T (8.47)

Second, we carry out the test by searching for the minimal t statistic for  = 1 hypothesis
in the following regression:
𝑘 𝑘 𝑘

ỹt = ∑ ω1i D(TB)1t−i + ∑ ω2i D(TB)2t−i + ỹt−1 + ∑ ci Δỹt−i + et

𝑖=0 𝑖=0 𝑖=1

t = k + 2, …., T (8.48)
where DUit = 1 if t > Tbi and 0 other wise, Tbi are the break dates, D(TB)it = 1 if t = Tbi
+ 1 and 0 otherwise.

86
Innovational outlier model:
The IO model is specified as:
Yt = μ + δ1 DU1t + δ2 DU2t + θ1 D(TB)1t + θ2 D(TB)2t
+Yt−1 + ∑𝑘𝑖=1 ci ΔYt−i + et (8.49)

8.5 Performing Unit Root Tests in Stata72

In this section, we use the example about the log(EX) series ([Link], Gujarati
(2011)) to illustrate the unit root tests using Stata. For all cases, the null hypothesis is
that log(EX) series represents a random walk, or has a unit root [except the case of
KPSS test as above mentioned]. We can learn the syntax in details and practice each of
the unit root tests online with webuse examples. In doing so, the commands are:
help dfuller
help pperron
help dfgls
help kpss
help zandrews
help clemao1
help clemio1
help clemao2
help clemio2
Table 8.4: DF test of log(EX).

The DF t statistic is 0.172, which is positive. This incorrect sign may be due to the test
equation is incorrectly specified. The positive  would imply  > 1, which means
log(EX) is explosive. This is not usual in macroeconomic data (Greene, 2008: p.740).
Therefore, we rule out this possibility.

72
For the first time, the commands such as kpss, dfgls, zandrews, clemao1, clemio1, clemeo2 or clemio2 are not
available. Therefore, you can use the command, say, findit kpss, then install it.
87
Table 8.5: DF test of log(EX).

The absolute value of the DF t statistic in this case is 3.026, less than the 10% critical
value of 3.128, so we should not reject the null hypothesis that log(EX) represents a
random walk, or has a unit root. In other words, log(EX) series is not stationary at 10%
level of significance. The MacKinnon approximate p-value of this test statistic is
approximately 12.48 percent as you can see at the bottom of the test results.

Table 8.6: DF-GLS test of log(EX).

The dfgls above reports tests of the nonstationary null hypothesis [i.e., the log(EX)
series represents a random walk, or has a unit root] for lags from 1 to 10 days. At the
bottom, the output offers three different methods for choosing an appropriate number
88
of lags: Ng-Perron sequential t, minimum Schwarz information criteria, and Ng-Perron
modified Akaike information criteria (MAIC). The MAIC is more recently developed,
and Monte Carlo experiments support its advantages over other criteria (Ng and Perron,
2001). The absolute value of DF-GLS statistic for 5 lags is 0.511, less than the 10%
critical value of 2.556, so we should not reject the null hypothesis. Note that Ng-Perron
sequential t indicates the maximum lags is 26. However, for space limited, we just
restrict at 10 maximum lag lengths.
Let use the maximum lags at 5 from DF-GLS test results, we find that the absolute value
of the ADF t statistic in Table 8.7 is 2.809, less than the 10% critical value of 3.120, so
we also fail reject the null hypothesis that log(EX) represents a random walk, or has a
unit root. As a result, we can conclude that log(EX) series is not stationary at 10% level
of significance.

Table 8.7: ADF test of log(EX).

Table 8.8: PP test of log(EX).

89
Table 8.9: KPSS test of log(EX) for trend stationarity hypothesis.

The absolute value of the PP t-statistic in Table 8.8 is 3.027, less than the 10% critical
value of 3.120, so we should not reject the null hypothesis that log(EX) represents a
random walk, or has a unit root. Therefore, both ADF and PP tests confirm that log(EX)
series is not stationary at 10% level of significance.
To make sure that the log(EX) series is not trend stationary, we use the KPSS test. The
results in Table 8.9 reject the null hypothesis that log(EX) is not trend stationary because
the test statistics at all lags are larger than critical values at 10%.

Table 8.10: dfgls test of log(EX).

A similar test of the first difference of log(EX) in Table 8.12, on the other hand, rejects
the nonstationary null hypothesis [i.e., log(EX) series has a unit root] at all lags [Note:
the maximum lag based on MAIC is up to 26], even at the 1% level. Therefore, we can
confirm that log(EX) series is a difference stationary process.
90
Table 8.11: ADF test of log(EX).

The absolute value of the ADF t-statistic is about 17.68 [Table 8.11] and PP t statistic
is about 48.39 [Table 8.12], greater than the 1% critical value of 3.43, so we should
reject the null hypothesis that the first difference of log(EX) series represents a random
walk, or has a unit root. Eventually, we could conclude that log(EX) series follows a
difference stationary process.
Table 8.12: PP test of log(EX).

Table 8.13: Zivot and Andrews test of log(EX) and log(EX).

log(EX) log(EX)
Model Conclusion
t-statistic Tb t-statistic Tb
A –3.60 957 –48.65*** 1263 I(1)
B –3.10 2002 –48.57*** 861 I(1)
C –3.70 822 –48.67*** 1263 I(1)
Notes. *** p < 1%. Stata commands: zandrews lnex, break(intercept) lagmethod(BIC), zandrews lnex,
break(trend) lagmethod(BIC), zandrews lnex, break(both) lagmethod(BIC); zandrews [Link],
break(intercept) lagmethod(BIC), zandrews [Link], break(trend) lagmethod(BIC), zandrews [Link],
break(both) lagmethod(BIC).
Table 8.13 indicates that the series ‘log(EX)’ and log(EX) exhibit an endogenous
structural trend, and the series ‘log(EX)’ is integrated of order 1. The estimated breaks

91
are different among models, so we expect that the extended Zivot and Andrews’s
framework could be applicable. This extended framework was proposed by Lumsdaine
and Papell (1997), which allow for two possible structural breaks (see Ling et al. (2013)
and Lumsdaine and Papell (1997)).
We use the webuse example to illustrate the Perron and Vogelsang (1992) and Clemente
et al. (1998) tests because the lnex series contains 2119 observations, so it takes long
time to obtain the test results.
webuse wpi1
clemao1 wpi
clemio1 wpi
clemio2 wpi, maxlag(6) graph
clemao1 [Link]
clemio1 [Link]
clemao2 [Link]
clemio2 [Link]
Table 8.14: P-V (1992) and Clemente et al. (1998) tests of the wpi series.
log(EX) log(EX)
Model Conclusion
t-statistic Tb t- Tb
statistic
AO1 –1.68 1981Q4 –2.86 1974Q1 I(?)
IO1 –1.86 1972Q3 –3.07 1974Q2 I(?)
AO2 –4.28 1975Q1,1980Q2 –2.24 1972Q2, 1980Q3 I(?)
IO2 –2.86 1972Q3,1978Q3 –6.56** 1972Q3, 1980Q4 I(1)
Notes. ** p < 5%.
Table 8.14 shows that the quarterly whole price index (pwi) seems to be stationary at
the first difference if we apply the Clemente et al. (1998) innovational outlier model.
The results from other test models are inconclusive. This implies that if the series have
structural breaks, we must be careful. Table 8.15 presents a typical Stata’s outcome of
the Clemente-Montanes-Reyes unit root test. Besides, we can add the option ‘graph’ in
the clemao or clemio command. The result is shown in Figure 8.16.
Table 8.15: Clemente et al. (1998) test with double mean shifts, IO model.

92
 Clemente-Montañés-Reyes double IO test for unit root
Test on wpi: breaks at 1972q3,1978q3
U.S. Wholesale Price Index
20 40 60 80 100 120

1960q1 1970q1 1980q1 1990q1

quarterly date

[Link]
4
2
[Link]
0
-2

1960q1 1970q1 1980q1 1990q1

quarterly date

Figure 8.16: Clemente et al. (1998) double IO test for unit root.

8.6 Unit Root Tests in Empirical Studies

The conventional tests of unit roots such as ADF, DF-GLS, PP, NP (Ng and Perron,
2001) and KPSS tests have been very popular in time series studies. Therefore, this
section just introduces some examples associated with tests of structural breaks.
Salim, Yao, Chen, and Zhang (2017) investigates the long-run relationship and short-
run dynamics between foreign investment (FDI) and energy consumption in China over
the 1982-2012 period. Four time series in their model include per capita energy
consumption (C), per capita real income (Y), trade openness (the sum of exports and
imports) (T), and the flow of FDI (F). For consistency in data definitions and collection
method, the authors use data from the World Development Indicator database.
Prior to cointegration and causality tests, they apply various tests for integration of the
variables such as ADF, DF-GLS and NP tests. However, the authors are aware of the
criticisms toward the conventional unit root tests due to their low size and power, as
well as the neglect of any structural breaks in the series. The structural breaks are based
on the past market reforms in China since 1978. As a result, they apply the Clemente-
Montanes-Reyes unit root test. The test results are shown in Table 8.16.
Table 8.16 presents the test results for either one or two structural breaks in the selected
variables. Both AO and IO models show that energy consumption, economic growth

93
and trade openness are I(1) variables, while foreign investment may be I(0). The results
also indicate that structural breaks in the level series consistently appear around the
years 1990 and 2001, which correspond to the accelerated reform post-Southern Tour
in 1992 and the accession of the World Trade Organization (WTO) in 2001.
Table 8.16: The Clemente-Montanes-Reyes structural break unit root test.

Tursoy and Faisal (2018) investigate the long-run interaction between stock prices, gold
prices and crude oil prices by applying monthly data from Turkey for the period between
January 1986 and November 2016. The data is obtained from the World Bank
Commodity Price Data (The Pink Sheet) for crude oil price (Brent, $/bbl) (CB) and gold
price ($/troy oz) (G), and the International Financial Statistics for stock prices (equity
price - EQ). Prior to implementing the cointegration test among the variables in the
model, they determine the order of integration by using the Perron and Vogelsang
(1992) and Clemente, Montanes and Reyes (1998) unit root tests because of structural
breaks. The results of the Perron and Vogelsang (1992) unit root test (Table 8.17) show
that all the variables in the series are nonstationary at level. However, the become
stationary by taking the first difference. The same inference is derived from the
Clemente, Montanes and Reyes (1998) unit root test (Table 8.18). Thus both unit root

94
tests demonstrate that the variables in the series have a unique order of integration,
which means that they are I(1).
Table 8.17: Perron-Vogelsang test with one endogenous structural break.

Table 8.18: Clemente-Montanes-Reyes test with two endogenous structural breaks.

Balaguer and Cantavella (2018) examines the role of education in the environmental
quality in Australia over the 1950-2014 period. Data of this study include per capita
fossil fuel emissions of CO2 (c), real income per capita (y), total number of higher
education students in the country at both graduate and postgraduate levels (d), and some
control variables such as real oil prices, Gini index, and trade openness. The authors
said that to overcome potential problems of identification under presence of structural
breaks, they use a modified Dickey Fuller unit root test (following Vogelsang and
Perron, 1998). The test results are presented in Table 8.19. The results from two
complementary tests (the additive outlier, AO, and innovational outlier, IO) are quite

95
consistent at the 5% level of significance. In particular, all series reject the presence of
two unit roots in favour of just one unit root (i.e., I(1) variables). In addition, the authors
conclude that the square of education time series in the IO test seems to be a trend
stationary process against a difference stationary process whereas trade would be an
I(0) variable. Therefore, a mix of I(0) and I(1) variables exists in the system.
Table 8.19: Modified augmented Dickey Fuller test (Vogelsang and Perron, 1998).

Shahbaz, Solarin, Hammoudeh, and Shahzad (2017) investigates the specification of

the environmental Kuznets curve for the US economy by accounting for the presence
of a major renewable energy source and trade openness over the period 1960-2016. Data
including real GDP per capita (Yt), real trade (exports + imports) per capita (TRt), CO2
emissions per capita (Ct) and oil prices (Ot) were collected from the World Development
Indicators, and biomass energy consumption per capita (Et) was collected from
[Link]. In order to examine the unit root properties of the variables, the
authors first use the NP test because it provides efficient empirical results for small
samples. The NP test results indicate that all the series are nonstationary in the level but
are stationary in the first difference of the variables. However, the authors said that the
Ng-Perron unit root test provides ambiguous empirical results due to their low
explanatory power since this unit root test does not accommodate information about
unknown structural break dates stemming from the series, which further weakens the
stationary hypothesis.
To resolve this issue, the authors employ the Clemente et al. (1998) unit root test which
contains information about single and double unknown structural breaks occurring in
the series during the sample period. The test results are presented in Table 8.20. The
results indicate that all the variable are nonstationary in the level, but stationary in their
first difference form in the presence of structural breaks. Therefore, they conclude that
all the variables are integrated of I(1).

96
Table 8.20: Clemente et al. (1998) unit root test results.

Fosten, Morley and Taylor (2012) investigates the behavior of emissions in

disequilibrium with respect to the environmental Kuznets curve relationship over the
1830-2003 period in the United Kingdom. The variables of interest include CO2 per
capita (et), SO2 per capita (et), real GDP per capita (yt), and energy price (ept). Before
performing the cointegration analysis, the authors check the order of integration of the
variables. First of all they run three basic tests with no structural breaks, including ADF,
DFGLS and PP. They then run the structural break unit root test of Zivot and Andrews
(1992) - ZA. The ZA test results are presented in Table 8.21.
Table 8.21: Zivot and Andrews (1992) test results.

97