This is page 57
Printer: Opaque this

3
Time Series Concepts

3.1 Introduction
This chapter provides background material on time series concepts that are used throughout the book. These concepts are presented in an informal way, and extensive examples using S-PLUS are used to build intuition. Section 3.2 discusses time series concepts for stationary and ergodic univariate time series. Topics include testing for white noise, linear and autoregressive moving average (ARMA) process, estimation and forecasting from ARMA models, and long-run variance estimation. Section 3.3 introduces univariate nonstationary time series and defines the important concepts of I(0) and
I(1) time series. Section 3.4 explains univariate long memory time series.
Section 3.5 covers concepts for stationary and ergodic multivariate time series, introduces the class of vector autoregression models, and discusses long-run variance estimation.
Rigorous treatments of the time series concepts presented in this chapter can be found in Fuller (1996) and Hamilton (1994). Applications of these concepts to financial time series are provided by Campbell, Lo and MacKinlay (1997), Mills (1999), Gourieroux and Jasiak (2001), Tsay (2001), Alexander (2001) and Chan (2002).

58

3. Time Series Concepts

3.2 Univariate Time Series
3.2.1 Stationary and Ergodic Time Series
Let {yt } = {. . . yt−1 , yt , yt+1 , . . .} denote a sequence of random variables indexed by some time subscript t. Call such a sequence of random variables a time series.
The time series {yt } is covariance stationary if
E[yt ] = µ for all t cov(yt , yt−j ) = E[(yt − µ)(yt−j − µ)] = γ j for all t and any j
For brevity, call a covariance stationary time series simply a stationary time series. Stationary time series have time invariant first and second moments. The parameter γ j is called the j th order or lag j autocovariance of {yt } and a plot of γ j against j is called the autocovariance function. The autocorrelations of {yt } are defined by γj cov(yt , yt−j ) ρj = p
=
γ0 var(yt )var(yt−j )

and a plot of ρj against j is called the autocorrelation function (ACF).
Intuitively, a stationary time series is defined by its mean, variance and
ACF. A useful result is that any function of a stationary time series is also a stationary time series. So if {yt } is stationary then {zt } = {g(yt )} is stationary for any function g(·).
The lag j sample autocovariance and lag j sample autocorrelation are defined as γj ˆ

=

ρj
ˆ

=

T
1 X
(yt − y )(yt−j − y )
¯
¯
T t=j+1

γj
ˆ
γ0
ˆ

(3.1)
(3.2)

PT
1
where y = T t=1 yt is the sample mean. The sample ACF (SACF) is a
¯
plot of ρj against j.
ˆ
A stationary time series {yt } is ergodic if sample moments converge in p p p ˆ
ˆ
probability to population moments; i.e. if y → µ, γ j → γ j and ρj → ρj .
¯
Example 1 Gaussian white noise (GWN) processes
Perhaps the most simple stationary time series is the independent Gaussian white noise process yt ∼ iid N (0, σ 2 ) ≡ GW N (0, σ 2 ). This process has µ = γ j = ρj = 0 (j 6= 0). To simulate a GW N (0, 1) process in S-PLUS use the rnorm function:

3.2 Univariate Time Series

59

-0.2

-2

0.0

-1

0.2

0

y

ACF
0.4

0.6

1

0.8

2

1.0

Series : y

0

20

40

60

80

100

0

2

4

Lag

6

8

10

FIGURE 3.1. Simulated Gaussian white noise process and SACF.

> set.seed(101)
> y = rnorm(100,sd=1)
To compute the sample moments y , γ j , ρj (j = 1, . . . , 10) and plot the
¯ ˆ ˆ data and SACF use
>
>
>
>
>
>

y.bar = mean(y)
g.hat = acf(y,lag.max=10,type="covariance",plot=F)
r.hat = acf(y,lag.max=10,type="correlation",plot=F) par(mfrow=c(1,2)) tsplot(y,ylab="y") acf.plot(r.hat) By default, as shown in Figure 3.1, the SACF is shown with 95% confidence limits about zero. These limits are based on the result (c.f. Fuller
(1996) pg. 336) that if {yt } ∼ iid (0, σ2 ) then µ ¶
1
A
, j > 0. ρj ∼ N 0,
ˆ
T

¡ 1¢
A
ˆ
The notation ρj ∼ N 0, T means that the distribution of ρj is approxiˆ
1
mated by normal distribution with mean 0 and variance T and is based on
√
d the central limit theorem result T ρj → N (0, 1). The 95% limits about
ˆ
1.96 zero are then ± √T .

3. Time Series Concepts

y

-2

-1

0

1

2

60

-2

-1

0

1

2

Quantiles of Standard Normal

FIGURE 3.2. Normal qq-plot for simulated GWN.

Two slightly more general processes are the independent white noise
(IWN) process, yt ∼ IW N (0, σ 2 ), and the white noise (WN) process, yt ∼ W N (0, σ 2 ). Both processes have mean zero and variance σ 2 , but the IWN process has independent increments, whereas the WN process has uncorrelated increments.
Testing for Normality
In the previous example, yt ∼ GW N (0, 1). There are several statistical methods that can be used to see if an iid process yt is Gaussian. The most common is the normal quantile-quantile plot or qq-plot, a scatterplot of the standardized empirical quantiles of yt against the quantiles of a standard normal random variable. If yt is normally distributed, then the quantiles will lie on a 45 degree line. A normal qq-plot with 45 degree line for yt may be computed using the S-PLUS functions qqnorm and qqline
> qqnorm(y)
> qqline(y)
Figure 3.2 shows the qq-plot for the simulated GWN data of the previous example. The quantiles lie roughly on a straight line. The S+FinMetrics function qqPlot may be used to create a Trellis graphics qq-plot.
The qq-plot is an informal graphical diagnostic. Two popular formal statistical tests for normality are the Shapiro-Wilks test and the Jarque-

3.2 Univariate Time Series

61

Bera test. The Shapiro-Wilk’s test is a well-known goodness of fit test for the normal distribution. It is attractive because it has a simple, graphical interpretation: one can think of it as an approximate measure of the correlation in a normal quantile-quantile plot of the data. The Jarque-Bera test is based on the result that a normally distributed random variable has skewness equal to zero and kurtosis equal to three. The Jarque-Bera test statistic is
Ã
! d T [ 2 (kurt − 3)2
JB = skew +
(3.3)
6
4

d
[
where skew denotes the sample skewness and kurt denotes the sample kurtosis. Under the null hypothesis that the data is normally distributed
A

JB ∼ χ2 (2).

Example 2 Testing for normality using the S+FinMetrics function normalTest The Shapiro-Wilks and Jarque-Bera statistics may be computed using the S+FinMetrics function normalTest. For the simulated GWN data of the previous example, these statistics are
> normalTest(y, method="sw")
Test for Normality: Shapiro-Wilks
Null Hypothesis: data is normally distributed
Test Statistics:
Test Stat 0.9703
p.value 0.1449
Dist. under Null: normal
Total Observ.: 100
> normalTest(y, method="jb")
Test for Normality: Jarque-Bera
Null Hypothesis: data is normally distributed
Test Statistics:
Test Stat 1.8763
p.value 0.3914

62

3. Time Series Concepts

Dist. under Null: chi-square with 2 degrees of freedom
Total Observ.: 100
The null of normality is not rejected using either test.
Testing for White Noise
Consider testing the null hypothesis
H0 : yt ∼ W N (0, σ 2 ) against the alternative that yt is not white noise. Under the null, all of the autocorrelations ρj for j > 0 are zero. To test this null, Box and Pierce
(1970) suggested the Q-statistic
Q(k) = T

k
X

ρ2
ˆj

(3.4)

j=1

where ρj is given by (3.2). Under the null, Q(k) is asymptotically disˆ tributed χ2 (k). In a finite sample, the Q-statistic (3.4) may not be well approximated by the χ2 (k). Ljung and Box (1978) suggested the modified
Q-statistic
k
X ρ2
ˆj
M Q(k) = T (T + 2)
(3.5)
T −j j=1 which is better approximated by the χ2 (k) in finite samples.
Example 3 Daily returns on Microsoft

Consider the time series behavior of daily continuously compounded returns on Microsoft for 2000. The following S-PLUS commands create the data and produce some diagnostic plots:
>
>
>
+
>
>
>
>
>

r.msft = getReturns(DowJones30[,"MSFT"],type="continuous")
r.msft@title = "Daily returns on Microsoft" sample.2000 = (positions(r.msft) > timeDate("12/31/1999")
& positions(r.msft) < timeDate("1/1/2001")) par(mfrow=c(2,2)) plot(r.msft[sample.2000],ylab="r.msft")
r.acf = acf(r.msft[sample.2000]) hist(seriesData(r.msft)) qqnorm(seriesData(r.msft))

The daily returns on Microsoft resemble a white noise process. The qqplot, however, suggests that the tails of the return distribution are fatter than the normal distribution. Notice that since the hist and qqnorm functions do not have methods for “timeSeries” objects the extractor function seriesData is required to extract the data frame from the data slot of r.msft.

3.2 Univariate Time Series
Series : r.msft[sample.2000]

-0.2

-0.15

0.2

0.05

ACF
0.4 0.6

0.8

0.15

1.0

Daily returns on Microsoft

r.msft

63

Dec
2001

0

seriesData(r.msft)

600
400

10

Lag

15

20

0

200

5

0.1

Oct

0.0

Jun
Aug
2000

-0.1

Apr

800

Feb

-0.1

0.0

0.1

seriesData(r.msft)

-2

0

2

Quantiles of Standard Normal

FIGURE 3.3. Daily returns on Microsoft with diagnostic plots.

The S+FinMetrics functions histPlot and qqPlot will produce a histogram and qq-plot for a “timeSeries” object using Trellis graphics. For example, > histPlot(r.msft,strip.text="MSFT monthly return")
> qqPlot(r.msft,strip.text="MSFT monthly return")
However, Trellis plots cannot be displayed in a multipanel plot created using par.
The S+FinMetrics function autocorTest may be used to compute the
Q-statistic and modified Q-statistic to test the null that the returns on
Microsoft follow a white noise process:
> autocorTest(r.msft, lag.n=10, method="lb")
Test for Autocorrelation: Ljung-Box
Null Hypothesis: no autocorrelation
Test Statistics:
Test Stat 11.7746
p.value 0.3004
Dist. under Null: chi-square with 10 degrees of freedom

64

3. Time Series Concepts

Total Observ.: 2527
The argument lag.n=10 specifies that k = 10 autocorrelations are used in computing the statistic, and method="lb" specifies that the modified
Box-Pierce statistic (3.5) be computed. To compute the simple Box-Pierce statistic, specify method="bp". The results indicate that the white noise null cannot be rejected.

3.2.2 Linear Processes and ARMA Models
Wold’s decomposition theorem (c.f. Fuller (1996) pg. 96) states that any covariance stationary time series {yt } has a linear process or infinite order moving average representation of the form yt = µ+

ψ0

= 1,

∞
X

ψ k εt−k

(3.6)

k=0
∞
X

k=0

ψ2 < ∞ k ∼ W N (0, σ 2 )

εt

In the Wold form, it can be shown that
E[yt ] = µ γ0 = var(yt ) = σ 2

∞
X

ψ2 k k=0

γj ρj = cov(yt , yt−j ) = σ2
=

P∞

k=0 ψ k ψ k+j
P∞
2 k=0 ψ k

∞
X

ψ k ψ k+j

k=0

Hence, the pattern of autocorrelations in any stationary and ergodic time series {yt } is determined by the moving average weights {ψ j } in its Wold representation. To ensure convergence of the linear process representation to a stationary and ergodic process with nice properties, it is necessary to further restrict the behavior of the moving average weights {ψ j }. A standard assumption used in the econometrics literature (c.f. Hamilton
(1994) pg. 504) is 1-summability
∞
X j=0 j|ψ j | = 1 + 2|ψ 2 | + 3|ψ 3 | + · · · < ∞.

The moving average weights in the Wold form are also called impulse responses since
∂yt+s
= ψ s , s = 1, 2, . . .
∂εt

3.2 Univariate Time Series

65

For a stationary and ergodicP time series lims→∞ ψ s = 0 and the long-run cumulative impulse response ∞ ψ s < ∞. A plot of ψ s against s is called s=0 the impulse response function (IRF).
The general Wold form of a stationary and ergodic time series is handy for theoretical analysis but is not practically useful for estimation purposes.
A very rich and practically useful class of stationary and ergodic processes is the autoregressive-moving average (ARMA) class of models made popular by Box and Jenkins (1976). ARMA(p, q) models take the form of a pth order stochastic difference equation yt − µ = φ1 (yt−1 − µ) + · · · + φp (yt−p − µ)
+εt + θ1 εt−1 + · · · + θq εt−q εt ∼ W N (0, σ 2 )

(3.7)

ARMA(p, q) models may be thought of as parsimonious approximations to the general Wold form of a stationary and ergodic time series. More information on the properties of ARMA(p, q) process and the procedures for estimating and forecasting these processes using S-PLUS are in the SPLUS Guide to Statistics Vol. II, chapter 27, Venables and Ripley (1999) chapter 13, and Meeker (2001)1 .
Lag Operator Notation
The presentation of time series models is simplified using lag operator notation. The lag operator L is defined such that for any time series {yt },
Lyt = yt−1 . The lag operator has the following properties: L2 yt = L · Lyt = yt−2 , L0 = 1 and L−1 yt = yt+1 . The operator ∆ = 1 − L creates the first difference of a time series: ∆yt = (1 − L)yt = yt − yt−1 . The ARMA(p, q) model (3.7) may be compactly expressed using lag polynomials. Define φ(L) = 1 − φ1 L − · · · − φp Lp and θ(L) = 1 + θ1 L + · · · + θq Lq . Then (3.7) may be expressed as φ(L)(yt − µ) = θ(L)εt
Similarly, the Wold representation in lag operator notation is yt = µ + ψ(L)εt
∞
X ψ(L) = ψ k Lk , ψ 0 = 1 k=0 and the long-run cumulative impulse response is ψ(1) (i.e. evaluate ψ(L) at L = 1). With ARMA(p, q) models the Wold polynomial ψ(L) is approx1 William Meeker also has a library of time series functions for the analysis of
ARMA models available for download at http://www.public.iastate.edu/~stat451/splusts/splusts.html. 66

3. Time Series Concepts

imated by the ratio of the AR and MA polynomials ψ(L) =

θ(L) φ(L) 3.2.3 Autoregressive Models
AR(1) Model
A commonly used stationary and ergodic time series in financial modeling is the AR(1) process yt − µ = φ(yt−1 − µ) + εt , t = 1, . . . , T where εt ∼ W N (0, σ 2 ) and |φ| < 1. The above representation is called the mean-adjusted form. The characteristic equation for the AR(1) is φ(z) = 1 − φz = 0

(3.8)

1 so that the root is z = φ . Stationarity is satisfied provided the absolute value of the root of the characteristic equation (3.8) is greater than one:
1
σ2
| φ | > 1 or |φ| < 1. In this case, it is easy to show that E[yt ] = µ, γ 0 = 1−φ2 ,

ψ j = ρj = φj and the Wold representation is yt = µ +

∞
X

ρj εt−j .

j=0

Notice that for the AR(1) the ACF and IRF are identical. This is not true
1
in general. The long-run cumulative impulse response is ψ(1) = 1−φ .
The AR(1) model may be re-written in components form as yt ut

= µ + ut
= φut−1 + εt

or in autoregression form as yt = c + φyt−1 + εt c = µ(1 − φ)
An AR(1) with µ = 1, φ = 0.75, σ 2 = 1 and T = 100 is easily simulated in S-PLUS using the components form:
>
>
>
>

set.seed(101) e = rnorm(100,sd=1)
e.start = rnorm(25,sd=1)
y.ar1 = 1 + arima.sim(model=list(ar=0.75), n=100,

3.2 Univariate Time Series
ACF and IRF for AR(1)

0.6
0.4

-2

0.2

0

2

Autocorrelation

4

Simulated AR(1)

67

0

20

40

60

80

100

1

2

3

4

5

6

7

8

9

10

lag

-0.2

0.2

ACF

0.6

1.0

Series : y.ar1

0

2

4

Lag

6

8

10

FIGURE 3.4. Simulated AR(1), ACF, IRF and SACF.

+ innov=e, start.innov=e.start)
> mean(y.ar1)
[1] 1.271
> var(y.ar1)
[1] 2.201
The ACF and IRF may be computed as
> gamma.j = rep(0.75,10)^seq(10)
The simulated data, ACF and SACF are illustrated in Figure 3.4 using
>
>
>
>
+
>

par(mfrow=c(2,2)) tsplot(y.ar1,main="Simulated AR(1)") abline(h=1) tsplot(gamma.j, type="h", main="ACF and IRF for AR(1)", ylab="Autocorrelation", xlab="lag") tmp = acf(y.ar1, lag.max=10)

Notice that {yt } exhibits mean-reverting behavior. That is, {yt } fluctuates about the mean value µ = 1. The ACF and IRF decay at a geometric rate.
The decay rate of the IRF is sometimes reported as a half-life — the lag j half at which the IRF reaches 1 . For the AR(1) with positive φ, it can be
2
shown that j half = ln(0.5)/ ln(φ). For φ = 0.75, the half-life is
> log(0.5)/log(0.75)

68

3. Time Series Concepts

-0.5

-0.1

0.4

US/CA 30 day interest rate differential

1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996

ACF

0.0

0.2

0.4

0.6

0.8

1.0

Series : uscn.id

0

5

10

15

20

Lag

FIGURE 3.5. US/CA 30 day interest rate differential and SACF.

[1] 2.409
Many economic and financial time series are well characterized by an
AR(1) process. Leading examples in finance are valuation ratios (dividendprice ratio, price-earning ratio etc), real exchange rates, interest rates, and interest rate differentials (spreads). To illustrate, consider the 30day US/CA interest rate differential2 constructed from the S+FinMetrics
“timeSeries” object lexrates.dat:
>
+
>
>
>
>
>
>

uscn.id = 100*(lexrates.dat[,"USCNF"]lexrates.dat[,"USCNS"]) colIds(uscn.id) = "USCNID" uscn.id@title = "US/CA 30 day interest rate differential" par(mfrow=c(2,1)) plot(uscn.id,reference.grid=F) abline(h=0) tmp = acf(uscn.id)

The interest rate differential is clearly persistent: autocorrelations are significant at the 5% level up to 15 months.
2 By covered interest rate parity, the nominal interest rate differential between risk free bonds from two countries is equal to the difference between the nominal forward and spot exchange rates.

3.2 Univariate Time Series

69

AR(p) Models
The AR(p) model in mean-adjusted form is yt − µ = φ1 (yt−1 − µ) + · · · + φp (yt−p − µ) + εt or, in lag operator notation, φ(L)(yt − µ) = εt where φ(L) = 1 − φ1 L − · · · − φp Lp . The autoregressive form is φ(L)yt = c + εt .
It can be shown that the AR(p) is stationary and ergodic provided the roots of the characteristic equation φ(z) = 1 − φ1 z − φ2 z 2 − · · · − φp z p = 0

(3.9)

lie outside the complex unit circle (have modulus greater than one). A necessary condition for stationarity that is useful in practice is that |φ1 +
· · · + φp | < 1. If (3.9) has complex roots then yt will exhibit sinusoidal behavior. In the stationary AR(p), the constant in the autoregressive form is equal to µ(1 − φ1 − · · · − φp ).
The moments of the AR(p) process satisfy the Yule-Walker equations γ0 γj

= φ1 γ 1 + φ2 γ 2 + · · · + φp γ p + σ 2
= φ1 γ j−1 + φ2 γ j−2 + · · · + φp γ j−p

(3.10)

A simple recursive algorithm for finding the Wold representation is based on matching coefficients in φ(L) and ψ(L) such that φ(L)ψ(L) = 1. For example, in the AR(2) model
(1 − φ1 L − φ2 L2 )(1 + ψ 1 L + ψ 2 L2 + · · · ) = 1 implies ψ1 ψ2 ψ3

ψj

= 1
= φ1 ψ 1 + φ2
= φ1 ψ 2 + φ2 ψ 1
.
.
.
= φ1 ψ j−1 + φ2 ψ j−2

Partial Autocorrelation Function
The partial autocorrelation function (PACF) is a useful tool to help identify AR(p) models. The PACF is based on estimating the sequence of AR

70

3. Time Series Concepts
Series : irate.real

0.0

-0.004

0.2

0.000

ACF
0.4 0.6

0.8

0.004

1.0

Monthly Real Interest Rate

0

5

10

Lag

15

20

25

1965

1970

1975

1980

1985

1990

1995

2000

0.0

Partial ACF
0.2 0.4 0.6

0.8

Series : irate.real

0

5

10

Lag

15

20

25

FIGURE 3.6. Monthly U.S. real interest rate, SACF and SPACF.

models zt zt

zt

= φ11 zt−1 + ε1t
= φ21 zt−1 + φ22 zt−2 + ε2t
.
.
.
= φp1 zt−1 + φp2 zt−2 + · · · + φpp zt−p + εpt

where zt = yt −µ is the demeaned data. The coefficients φjj for j = 1, . . . , p
(i.e., the last coefficients in each AR(p) model) are called the partial autocorrelation coefficients. In an AR(1) model the first partial autocorrelation coefficient φ11 is non-zero, and the remaining partial autocorrelation coefficients φjj for j > 1 are equal to zero. Similarly, in an AR(2), the first and second partial autocorrelation coefficients φ11 and φ22 are non-zero and the rest are zero for j > 2. For an AR(p) all of the first p partial autocorrelation coefficients are non-zero, and the rest are zero for j > p.
The sample partial autocorrelation coefficients up to lag p are essentially obtained by estimating the above sequence of p AR models by least squares
ˆ
and retaining the estimated coefficients φjj .
Example 4 Monthly real interest rates

3.2 Univariate Time Series

71

The “timeSeries” object varex.ts in the S+FinMetrics module contains monthly data on real stock returns, real interest rates, inflation and real output growth.
> colIds(varex.ts)
[1] "MARKET.REAL" "RF.REAL"

"INF"

"IPG"

Figure 3.6 shows the real interest rate, RF.REAL, over the period January
1961 through December 2000 produced with the S-PLUS commands
>
>
>
>
>
>

smpl = (positions(varex.ts) > timeDate("12/31/1960")) irate.real = varex.ts[smpl,"RF.REAL"] par(mfrow=c(2,2)) acf.plot(acf(irate.real, plot=F)) plot(irate.real, main="Monthly Real Interest Rate") tmp = acf(irate.real, type="partial")

The SACF and SPACF indicate that the real interest rate might be modeled as an AR(2) or AR(3) process.

3.2.4 Moving Average Models
MA(1) Model
The MA(1) model has the form yt = µ + εt + θεt−1 , εt ∼ W N (0, σ 2 )
For any finite θ the MA(1) is stationary and ergodic. The moments are
E[yt ] = µ, γ 0 = σ 2 (1 + θ2 ), γ 1 = σ2 θ, γ j = 0 for j > 1 and ρ1 = θ/(1 + θ2 ).
Hence, the ACF of an MA(1) process cuts off at lag one, and the maximum value of this correlation is ±0.5.
There is an identification problem with the MA(1) model since θ = 1/θ produce the same value of ρ1 . The MA(1) is called invertible if |θ| < 1 and is called non-invertible if |θ| ≥ 1. In the invertible MA(1), the error term εt has an infinite order AR representation of the form εt =

∞
X
j=0

θ∗j (yt−j − µ)

where θ∗ = −θ so that εt may be thought of as a prediction error based on past values of yt . A consequence of the above result is that the PACF for an invertible MA(1) process decays towards zero at an exponential rate.
Example 5 Signal plus noise model

72

3. Time Series Concepts
1st difference

0

-3

-2

-2

-1

y

dy

0

2

1

2

Signal plus noise

0

20

40

60

80

100

0

20

40

80

100

Series : dy

-0.5

-0.4

0.0

ACF

0.5

Partial ACF
-0.2
0.0

1.0

0.2

Series : dy

60

0

5

10
Lag

15

0

5

10
Lag

15

FIGURE 3.7. Simulated data, SACF and SPACF from signal plus noise model.

MA(1) models often arise through data transformations like aggregation and differencing3 . For example, consider the signal plus noise model yt zt

= zt + εt , εt ∼ W N (0, σ 2 ) ε = zt−1 + ηt , η t ∼ W N (0, σ2 ) η where εt and η t are independent. For example, zt could represent the fundamental value of an asset price and εt could represent an iid deviation about the fundamental price. A stationary representation requires differencing yt :
∆yt = ηt + εt − εt−1
It can be shown, e.g. Harvey (1993), that ∆yt is an MA(1) process with θ =
√
−(q+2)+
2

q 2 +4q

σ2

−1 where q = ση is the signal-to-noise ratio and ρ1 = q+2 < 0.
2
ε
Simulated data with σ 2 = 1 and σ 2 = (0.5)2 created with the S-PLUS ε η commands > set.seed(112)
> eps = rnorm(100,sd=1)
> eta = rnorm(100,sd=0.5)
3 MA(1) type models for asset returns often occur as the result of no-trading effects or bid-ask bounce effects. See Campbell, Lo and MacKinlay (1997) chapter 3 for details.

3.2 Univariate Time Series

>
>
>
>
>
>
>
>

73

z = cumsum(eta) y = z + eps dy = diff(y) par(mfrow=c(2,2)) tsplot(y, main="Signal plus noise",ylab="y") tsplot(dy, main="1st difference",ylab="dy") tmp = acf(dy) tmp = acf(dy,type="partial")

are illustrated in Figure 3.7. The signal-to-noise ratio q = 0.25 implies a first lag autocorrelation of ρ1 = −0.444. This negative correlation is clearly reflected in the SACF.
MA(q) Model
The MA(q) model has the form yt = µ + εt + θ1 εt−1 + · · · + θq εt−q , where εt ∼ W N (0, σ 2 )
The MA(q) model is stationary and ergodic provided θ1 , . . . , θ q are finite.
It is invertible if all of the roots of the MA characteristic polynomial θ(z) = 1 + θ1 z + · · · θq z q = 0

(3.11)

lie outside the complex unit circle. The moments of the MA(q) are
E[yt ] = µ γ 0 = σ 2 (1 + θ2 + · · · + θ2 )
1
q
½
(θj + θj+1 θ1 + θj+2 θ2 + · · · + θq θq−j ) σ2 for j = 1, 2, . . . , q γj =
0 for j > q
Hence, the ACF of an MA(q) is non-zero up to lag q and is zero afterwards.
As with the MA(1), the PACF for an invertible MA(q) will show exponential decay and possibly pseudo cyclical behavior if the roots of (3.11) are complex. Example 6 Overlapping returns and MA(q) models
MA(q) models often arise in finance through data aggregation transformations. For example, let Rt = ln(Pt /Pt−1 ) denote the monthly continuously compounded return on an asset with price Pt . Define the annual return at time t using monthly returns as Rt (12) = ln(Pt /Pt−12 ) =
P11
2 j=0 Rt−j . Suppose Rt ∼ W N (µ, σ ) and consider a sample of monthly returns of size T , {R1 , R2 , . . . , RT }. A sample of annual returns may be created using overlapping or non-overlapping returns. Let {R12 (12), R13 (12),
. . . , RT (12)} denote a sample of T ∗ = T − 11 monthly overlapping annual returns and {R12 (12), R24 (12), . . . , RT (12)} denote a sample of T /12 nonoverlapping annual returns. Researchers often use overlapping returns in

74

3. Time Series Concepts

analysis due to the apparent larger sample size. One must be careful using overlapping returns because the monthly annual return sequence {Rt (12)} is not a white noise process even if the monthly return sequence {Rt } is.
To see this, straightforward calculations give
E[Rt (12)] γ0 γj γj =
=
=
=

12µ var(Rt (12)) = 12σ 2 cov(Rt (12), Rt−j (12)) = (12 − j)σ 2 for j < 12
0 for j ≥ 12

Since γ j = 0 for j ≥ 12 notice that {Rt (12)} behaves like an MA(11) process Rt (12) = 12µ + εt + θ1 εt−1 + · · · + θ11 εt−11 εt ∼ W N (0, σ 2 )
To illustrate, consider creating annual overlapping continuously compounded returns on the S&P 500 index over the period February 1990 through January 2001. The S+FinMetrics “timeSeries” singleIndex.dat contains the S&P 500 price data and the continuously compounded monthly returns are computed using the S+FinMetrics function getReturns
> sp500.mret = getReturns(singleIndex.dat[,"SP500"],
+ type="continuous")
> sp500.mret@title = "Monthly returns on S&P 500 Index"
The monthly overlapping annual returns are easily computed using the
S-PLUS function aggregateSeries
> sp500.aret = aggregateSeries(sp500.mret,moving=12,FUN=sum)
> sp500.aret@title = "Monthly Annual returns on S&P 500 Index"
The optional argument moving=12 specifies that the sum function is to be applied to moving blocks of size 12. The data together with the SACF and SPACF of the monthly annual returns are displayed in Figure 3.8.
The SACF has non-zero values up to lag 11. Interestingly, the SPACF is very small at all lags except the first.

3.2.5 ARMA(p,q) Models
The general ARMA(p, q) model in mean-adjusted form is given by (3.7).
The regression formulation is yt = c + φ1 yt−1 + · · · + φp yt−p + εt + θεt−1 + · · · + θεt−q

(3.12)

It is stationary and ergodic if the roots of the characteristic equation φ(z) =
0 lie outside the complex unit circle, and it is invertible if the roots of the

3.2 Univariate Time Series

Monthly Annual returns on S&P 500 Index

-0.10

-0.15

0.10

-0.05

0.05

0.30

Monthly returns on S&P 500 Index

75

1990

1992

1994

1996

1998

2000

1990

1994

1996

1998

2000

Series : sp500.aret

-0.2

-0.2

0.0

0.2

Partial ACF
0.2 0.4 0.6

ACF
0.4 0.6 0.8

0.8

1.0

Series : sp500.aret

1992

0

5

10
Lag

15

20

0

5

10
Lag

15

20

FIGURE 3.8. Monthly non-overlapping and overlapping annual returns on the
S&P 500 index.

MA characteristic polynomial θ(z) = 0 lie outside the unit circle. It is assumed that the polynomials φ(z) = 0 and θ(z) = 0 do not have canceling or common factors. A stationary and ergodic ARMA(p, q) process has a mean equal to c µ=
(3.13)
1 − φ1 − · · · − φp and its autocovariances, autocorrelations and impulse response weights satisfy the recursive relationships γj ρj ψj = φ1 γ j−1 + φ2 γ j−2 + · · · + φp γ j−p
= φ1 ρj−1 + φ2 ρj−2 + · · · + φp ρj−p
= φ1 ψ j−1 + φ2 ψ j−2 + · · · + φp ψ j−p

The general form of the ACF for an ARMA(p, q) process is complicated.
See Hamilton (1994) chapter five for details. In general, for an ARMA(p, q) process, the ACF behaves like the ACF for an AR(p) process for p > q, and the PACF behaves like the PACF for an MA(q) process for q > p. Hence, both the ACF and PACF eventually show exponential decay.
ARMA(p, q) models often arise from certain aggregation transformations of simple time series models. An important result due to Granger and Morris (1976) is that if y1t is an ARMA(p1 , q1 ) process and y2t is an ARMA(p2 , q2 ) process, which may be contemporaneously correlated

76

3. Time Series Concepts

with y1t , then y1t + y2t is an ARMA(p, q) process with p = p1 + p2 and q = max(p1 + q2 , q1 + p2 ). For example, if y1t is an AR(1) process and y2 is a AR(1) process, then y1 + y2 is an ARMA(2,1) process.
High order ARMA(p, q) processes are difficult to identify and estimate in practice and are rarely used in the analysis of financial data. Low order
ARMA(p, q) models with p and q less than three are generally sufficient for the analysis of financial data.
ARIMA(p, d, q) Models
The specification of the ARMA(p, q) model (3.7) assumes that yt is stationary and ergodic. If yt is a trending variable like an asset price or a macroeconomic aggregate like real GDP, then yt must be transformed to stationary form by eliminating the trend. Box and Jenkins (1976) advocate removal of trends by differencing. Let ∆ = 1−L denote the difference operator. If there is a linear trend in yt then the first difference ∆yt = yt − yt−1 will not have a trend. If there is a quadratic trend in yt , then ∆yt will contain a linear trend but the second difference ∆2 yt = (1 − 2L + L2 )yt = yt − 2yt−1 + yt−2 will not have a trend. The class of ARMA(p, q) models where the trends have been transformed by differencing d times is denoted ARIMA(p, d, q)4 .

3.2.6 Estimation of ARMA Models and Forecasting
ARMA(p, q) models are generally estimated using the technique of maximum likelihood, which is usually accomplished by putting the ARMA(p, q) in state-space form from which the prediction error decomposition of the log-likelihood function may be constructed. Details of this process are given in Harvey (1993). An often ignored aspect of the maximum likelihood estimation of ARMA(p, q) models is the treatment of initial values. These initial values are the first p values of yt and q values of εt in (3.7). The exact likelihood utilizes the stationary distribution of the initial values in the construction of the likelihood. The conditional likelihood treats the p initial values of yt as fixed and often sets the q initial values of εt to zero. The exact maximum likelihood estimates (MLEs) maximize the exact log-likelihood, and the conditional MLEs maximize the conditional log-likelihood. The exact and conditional MLEs are asymptotically equivalent but can differ substantially in small samples, especially for models that are close to being nonstationary or noninvertible.5
4 More general ARIMA(p, d, q) models allowing for seasonality are discussed in chapter
27 of the S-PLUS Guide to Statistics, Vol. II.
5 As pointed out by Venables and Ripley (1999) page 415, the maximum likelihood estimates computed using the S-PLUS function arima.mle are conditional MLEs. Exact
MLEs may be easily computed using the S+FinMetrics state space modeling functions.

3.2 Univariate Time Series

77

For pure AR models, the conditional MLEs are equivalent to the least squares estimates from the model yt = c + φ1 yt−1 + · · · + φp yt−p + εt

(3.14)

Notice, however, that c in (3.14) is not an estimate of E[yt ] = µ. The least squares estimate of µ is given by plugging in the least squares estimates of c, φ1 , . . . , φp into (3.13).
Model Selection Criteria
Before an ARMA(p, q) may be estimated for a time series yt , the AR and
MA orders p and q must be determined by visually inspecting the SACF and SPACF for yt . Alternatively, statistical model selection criteria may be used. The idea is to fit all ARMA(p, q) models with orders p ≤ pmax and q ≤ qmax and choose the values of p and q which minimizes some model selection criteria. Model selection criteria for ARMA(p, q) models have the form M SC(p, q) = ln(˜ 2 (p, q)) + cT · ϕ(p, q) σ where σ 2 (p, q) is the MLE of var(εt ) = σ 2 without a degrees of freedom cor˜ rection from the ARMA(p, q) model, cT is a sequence indexed by the sample size T , and ϕ(p, q) is a penalty function which penalizes large ARMA(p, q) models. The two most common information criteria are the Akaike (AIC) and Schwarz-Bayesian (BIC):
2
(p + q)
T
ln T
BIC(p, q) = ln(˜ 2 (p, q)) + σ (p + q)
T
AIC(p, q) = ln(˜ 2 (p, q)) + σ The AIC criterion asymptotically overestimates the order with positive probability, whereas the BIC estimate the order consistently under fairly general conditions if the true orders p and q are less than or equal to pmax and qmax . However, in finite samples the BIC generally shares no particular advantage over the AIC.
Forecasting Algorithm
Forecasts from an ARIMA(p, d, q) model are straightforward. The model is put in state space form, and optimal h-step ahead forecasts along with forecast standard errors (not adjusted for parameter uncertainty) are produced using the Kalman filter algorithm. Details of the method are given in Harvey (1993).

78

3. Time Series Concepts

Estimation and Forecasting ARIMA(p, d, q) Models Using the S-PLUS
Function arima.mle
Conditional MLEs may be computed using the S-PLUS function arima.mle.
The form of the ARIMA(p, d, q) assumed by arima.mle is yt = φ1 yt−1 + · · · + φp yt−p
+εt − θ1 εt−1 − · · · − θq εt−q
+β0 xt

where xt represents additional explanatory variables. It is assumed that yt has been differenced d times to remove any trends and that the unconditional mean µ has been subtracted out so that yt is demeaned. Notice that arima.mle assumes that the signs on the MA coefficients θj are the opposite to those in (3.7).
The arguments expected by arima.mle are
> args(arima.mle) function(x, model = NULL, n.cond = 0, xreg = NULL, ...) where x is a univariate “timeSeries” or vector, model is a list object describing the specification of the ARMA model, n.cond sets the number of initial observations on which to condition in the formation of the log-likelihood, and xreg is a “timeSeries”, vector or matrix of additional explanatory variables. By default, arima.mle assumes that the
ARIMA(p, d, q) model is stationary and in mean-adjusted form with an estimate of µ subtracted from the observed data yt . To estimate the regression form (3.12) of the ARIMA(p, q) model, simply set xreg=1. ARIMA(p, d, q) models are specified using list variables the form
> mod.list = list(order=c(1,0,1))
> mod.list = list(order=c(1,0,1),ar=0.75,ma=0)
> mod.list = list(ar=c(0.75,-0.25),ma=c(0,0))
The first list simply specifies an ARMA(1,0,1)/ARMA(1,1) model. The second list specifies an ARIMA(1,0,1) as well as starting values for the
AR and MA parameters φ and θ. The third list implicitly determines an
ARMA(2,2) model by giving the starting values for the AR and MA parameters. The function arima.mle produces an object of class “arima” for which there are print and plot methods. Diagnostics from the fit can be created with the S-PLUS function arima.diag, and forecasts may be produced using arima.forecast.
Example 7 Estimation of ARMA model for US/CA interest rate differential
Consider estimating an ARMA(p, q) for the monthly US/CA interest rate differential data in the “timeSeries” uscn.id used in a previous

3.2 Univariate Time Series

79

example. To estimate an ARMA(1,1) model for the demeaned interest rate differential with starting values φ = 0.75 and θ = 0 use
> uscn.id.dm = uscn.id - mean(uscn.id)
> arma11.mod = list(ar=0.75,ma=0)
> arma11.fit = arima.mle(uscn.id.dm,model=arma11.mod)
> class(arma11.fit)
[1] "arima"
The components of arma11.fit are
> names(arma11.fit)
[1] "model"
"var.coef" "method"
"series"
[5] "aic"
"loglik"
"sigma2"
"n.used"
[9] "n.cond"
"converged" "conv.type" "call"
To see the basic fit simply type
> arma11.fit
Call: arima.mle(x = uscn.id.dm, model = arma11.mod)
Method: Maximum Likelihood
Model : 1 0 1
Coefficients:
AR : 0.82913
MA : 0.11008
Variance-Covariance Matrix: ar(1) ma(1) ar(1) 0.002046 0.002224 ma(1) 0.002224 0.006467
Optimizer has converged
Convergence Type: relative function convergence
AIC: -476.25563
ˆ
The conditional MLEs are φcmle = 0.829 and ˆcmle = −0.110. Standard θ errors for these parameters are given by the square roots of the diagonal elements of variance-covariance matrix
> std.errs = sqrt(diag(arma11.fit$var.coef))
> names(std.errs) = colIds(arma11.fit$var.coef)
> std.errs ar(1) ma(1)
0.04523 0.08041
It appears that the ˆcmle is not statistically different from zero. θ To estimate the ARMA(1,1) for the interest rate differential data in regression form (3.12) with an intercept use

80

3. Time Series Concepts

> arma11.fit2 = arima.mle(uscn.id,model=arma11.mod,xreg=1)
> arma11.fit2
Call: arima.mle(x = uscn.id, model = arma11.mod, xreg = 1)
Method: Maximum Likelihood
Model : 1 0 1
Coefficients:
AR : 0.82934
MA : 0.11065
Variance-Covariance Matrix: ar(1) ma(1) ar(1) 0.002043 0.002222 ma(1) 0.002222 0.006465
Coeffficients for regressor(s): intercept
[1] -0.1347
Optimizer has converged
Convergence Type: relative function convergence
AIC: -474.30852
The conditional MLEs for φ and θ are essentially the same as before, and the MLE for c is ccmle = −0.1347. Notice that the reported varianceˆ covariance matrix only gives values for the estimated ARMA coefficients
ˆ
φcmle and ˆcmle . θ Graphical diagnostics of the fit produced using the plot method
> plot(arma11.fit) are illustrated in Figure 3.9. There appears to be some high order serial correlation in the errors as well as heteroskedasticity.
The h-step ahead forecasts of future values may be produced with the
S-PLUS function arima.forecast. For example, to produce monthly forecasts for the demeaned interest rate differential from July 1996 through
June 1997 use
> fcst.dates = timeSeq("7/1/1996", "6/1/1997",
+ by="months", format="%b %Y")
> uscn.id.dm.fcst = arima.forecast(uscn.id.dm, n=12,
+ model=arma11.fit$model, future.positions=fcst.dates)
> names(uscn.id.dm.fcst)
[1] "mean"
"std.err"
The object uscn.id.dm.fcst is a list whose first component is a
“timeSeries” containing the h-step forecasts, and the second component is a “timeSeries” containing the forecast standard errors:
> uscn.id.dm.fcst[[1]]

3.2 Univariate Time Series

81

ARIMA Model Diagnostics: uscn.id.dm

-5

-1

3

Plot of Standardized Residuals

1978

1980

1982

1984 Plot of Residuals
1988
ACF 1986

1990

1992

1994

1996

-1.0

ACF
0.0

0.5

1.0

1976

0

5

10

15

20

-0.1

PACF
0.1

0.2

PACF Plot of Residuals

5

10

15

20

0.0

p-value
0.2
0.4

0.6

P-values of Ljung-Box Chi-Squared Statistics

4

6

8

10

12

14

Lag

ARIMA(1,0,1) Model with Mean 0

FIGURE 3.9. Residual diagnostics from ARMA(1,1) fit to US/CA interest rate differentials. Positions
Jul 1996
Aug 1996
Sep 1996
Oct 1996
Nov 1996
Dec 1996
Jan 1997
Feb 1997
Mar 1997
Apr 1997
May 1997
Jun 1997

1
0.09973
0.08269
0.06856
0.05684
0.04713
0.03908
0.03240
0.02686
0.02227
0.01847
0.01531
0.01270

The data, forecasts and 95% forecast confidence intervals shown in Figure
3.10 are produced by
>
>
+
+
+

smpl = positions(uscn.id.dm) >= timeDate("6/1/1995") plot(uscn.id.dm[smpl,],uscn.id.dm.fcst$mean, uscn.id.dm.fcst$mean+2*uscn.id.dm.fcst$std.err, uscn.id.dm.fcst$mean-2*uscn.id.dm.fcst$std.err, plot.args=list(lty=c(1,4,3,3)))

82

3. Time Series Concepts

-0.2

-0.1

0.0

0.1

0.2

0.3

US/CA 30 day interest rate differential

Q3

Q4

Q1

Q2

1995

Q3

Q4

1996

Q1

Q2
1997

FIGURE 3.10. Forecasts for 12 months for the series uscn.id.dm.

Estimating AR(p) by Least Squares Using the S+FinMetrics Function
OLS
As previously mentioned, the conditional MLEs for an AR(p) model may be computed using least squares. The S+FinMetrics function OLS, which extends the S-PLUS function lm to handle general time series regression, may be used to estimate an AR(p) in a particularly convenient way. The general use of OLS is discussed in Chapter 6, and its use for estimating an
AR(p) is only mentioned here. For example, to estimate an AR(2) model for the US/CA interest rate differential use
> ar2.fit = OLS(USCNID~ar(2), data=uscn.id)
> ar2.fit
Call:
OLS(formula = USCNID ~ar(2), data = uscn.id)
Coefficients:
(Intercept) lag1 -0.0265
0.7259

lag2
0.0758

Degrees of freedom: 243 total; 240 residual
Time period: from Apr 1976 to Jun 1996
Residual standard error: 0.09105

3.2 Univariate Time Series

83

ˆ
The least squares estimates of the AR coefficients are φ1 = 0.7259 and
ˆ
ˆ
ˆ
φ2 = 0.0758. Since φ1 + φ2 < 1 the estimated AR(2) model is stationary.
ˆ
ˆ
To be sure, the roots of φ(z) = 1 − φ1 z − φ2 z 2 = 0 are
> abs(polyroot(c(1,-ar2.fit$coef[2:3])))
[1] 1.222 10.798 are outside the complex unit circle.

3.2.7 Martingales and Martingale Difference Sequences
Let {yt } denote a sequence of random variables and let It = {yt, yt−1 , . . .} denote a set of conditioning information or information set based on the past history of yt . The sequence {yt , It } is called a martingale if
• It−1 ⊂ It (It is a filtration)
• E[|yt |] < ∞
• E[yt |It−1 ] = yt−1 (martingale property)
The most common example of a martingale is the random walk model yt = yt−1 + εt , εt ∼ W N (0, σ 2 ) where y0 is a fixed initial value. Letting It = {yt , . . . , y0 } implies E[yt |It−1 ] = yt−1 since E[εt |It−1 ] = 0.
Let {εt } be a sequence of random variables with an associated information set It . The sequence {εt , It } is called a martingale difference sequence
(MDS) if
• It−1 ⊂ It
• E[εt |It−1 ] = 0 (MDS property)
If {yt , It } is a martingale, a MDS {εt , It } may be constructed by defining εt = yt − E[yt |It−1 ]
By construction, a MDS is an uncorrelated process. This follows from the law of iterated expectations. To see this, for any k > 0
E[εt εt−k ] = E[E[εt εt−k |It−1 ]]
= E[εt−k E[εt |It−1 ]]
= 0
In fact, if zn is any function of the past history of εt so that zn ∈ It−1 then
E[εt zn ] = 0

84

3. Time Series Concepts

Although a MDS is an uncorrelated process, it does not have to be an independent process. That is, there can be dependencies in the higher order moments of εt . The autoregressive conditional heteroskedasticity (ARCH) process in the following example is a leading example in finance.
MDSs are particularly nice to work with because there are many useful convergence results (laws of large numbers, central limit theorems etc.).
White (1984), Hamilton (1994) and Hayashi (2000) describe the most useful of these results for the analysis of financial time series.
Example 8 ARCH process
A well known stylized fact about high frequency financial asset returns is that volatility appears to be autocorrelated. A simple model to capture such volatility autocorrelation is the ARCH process due to Engle (1982).
To illustrate, let rt denote the daily return on an asset and assume that
E[rt ] = 0. An ARCH(1) model for rt is rt zt σ2 t

= σ t zt
∼ iid N (0, 1)
2
= ω + αrt−1

(3.15)
(3.16)

where ω > 0 and 0 < α < 1. Let It = {rt , . . .}. The S+FinMetrics function simulate.garch may be used to generate simulations from above ARCH(1) model. For example, to simulate 250 observations on rt with ω = 0.1 and α = 0.8 use
> rt = simulate.garch(model=list(a.value=0.1, arch=0.8),
+
n=250, rseed=196)
> class(rt)
[1] "structure"
> names(rt)
[1] "et"
"sigma.t"
Notice that the function simulate.garch produces simulated values of both rt and σ t . These values are shown in Figure 3.11.
To see that {rt , It } is a MDS, note that
E[rt |It−1 ] = E[zt σ t |It−1 ]
= σ t E[zt |It−1 ]
= 0
Since rt is a MDS, it is an uncorrelated process. Provided |α| < 1, rt is a mean zero covariance stationary process. The unconditional variance of rt is given by
2
2 var(rt ) = E[rt ] = E[E[zt σ 2 |It−1 ]] t 2
2
= E[σ t E[zt |It−1 ] = E[σ 2 ] t 3.2 Univariate Time Series

85

-2

-1

0

1

2

Simulated returns

0

50

100

150

200

250

200

250

0.5

1.0

1.5

Simulated volatility

0

50

100

150

FIGURE 3.11. Simulated values from ARCH(1) process with ω = 1 and α = 0.8.
2
since E[zt |It−1 ] = 1. Utilizing (3.16) and the stationarity of rt , E[σ 2 ] may t be expressed as ω E[σ 2 ] = t 1−α

Furthermore, by adding ε2 to both sides of (3.16) and rearranging it follows t 2 that rt has an AR(1) representation of the form ε2 = ω + αε2 + vt t t−1

where vt = ε2 − σ 2 is a MDS. t t

3.2.8 Long-run Variance
Let yt be a stationary and ergodic time series. Anderson’s central limit theorem for stationary and ergodic processes (c.f. Hamilton (1994) pg. 195) states ∞
X
√ d T (¯ − µ) → N (0, y γj ) j=−∞ or




∞
1 X γ  y ∼ N µ,
¯
T j=−∞ j
A

86

3. Time Series Concepts

The sample size, T , times the asymptotic variance of the sample mean is often called the long-run variance of yt 6 : lrv(yt ) = T · avar(¯) = y ∞
X

γj .

j=−∞

Since γ −j = γ j , lrv(yt ) may be alternatively expressed as lrv(yt ) = γ 0 + 2

∞
X

γj .

j=1

Using the long-run variance, an asymptotic 95% confidence interval for µ takes the form q c y ± 1.96 · T −1 lrv(yt )
¯

c where lrv(yt ) is a consistent estimate of lrv(yt ).
Estimating the Long-Run Variance

If yt is a linear process, it may be shown that

2
∞
∞
X
X γ j = σ2  ψ j  = σ 2 ψ(1)2 j=−∞ j=0

and so

lrv(yt ) = σ 2 ψ(1)2

(3.17)

Further, if yt ∼ ARMA(p, q) then ψ(1) =

θ(1)
1 + θ1 + · · · + θq
=
1 − φ1 − · · · − φp φ(1) so that lrv(yt ) =

σ2 θ(1)2
.
φ(1)2

(3.18)

A consistent estimate of lrv(yt ) may then be computed by estimating the parameters of the appropriate ARMA(p, q) model and substituting these estimates into (3.18). Alternatively, the ARMA(p, q) process may be approximated by a high order AR(p∗ ) process yt = c + φ1 yt−1 + · · · + φp∗ yt−p∗ + εt
6 Using

spectral methods, lrv(¯) has the alternative representation y 1
2πf (0)
T
where f (0) denotes the spectral density of yt evaluated at frequency 0. lrv(¯) = y 3.2 Univariate Time Series

87

where the lag length p∗ is chosen such that εt is uncorrelated. This gives rise to the autoregressive long-run variance estimate lrvAR (yt ) =

σ2
.
φ (1)2
∗

(3.19)

A consistent estimate of lrv(yt ) may also be computed using some nonparametric methods. An estimator made popular by Newey and West
(1987) is the weighted autocovariance estimator c ˆ lrv NW (yt ) = γ 0 + 2

MT
X
j=1

wj,T · γ j
ˆ

(3.20)

where wj,T are weights which sum to unity and MT is a truncation lag parameter that satisfies MT = O(T 1/3 ). For MA(q) processes, γ j = 0 for j > q and Newey and West suggest using the rectangular weights wj,T = 1 for j ≤ MT = q; 0 otherwise. For general linear processes, Newey and West j suggest using the Bartlett weights wj,T = 1 − MT +1 with MT equal to the
2/9
integer part of 4(T /100) .
Example 9 Long-run variance of AR(1)
Let yt be an AR(1) process created using
> set.seed(101)
> e = rnorm(100,sd=1)
> y.ar1 = 1 + arima.sim(model=list(ar=0.75),innov=e)
Here ψ(1) =

1 φ(1) =

1
1−φ

and lrv(yt ) =

σ2
.
(1 − φ)2

For φ = 0.75, σ 2 = 1, lrv(yt ) = 16 implies for T = 100 an asymptotic standard error for y equal to SE(¯) = 0.40. If yt ∼ W N (0, 1), then the
¯
y asymptotic standard error for y is SE(¯) = 0.10.
¯
y lrvAR (yt ) may be easily computed in S-PLUS using OLS to estimate the
AR(1) parameters:
> ar1.fit = OLS(y.ar1~ar(1))
> rho.hat = coef(ar1.fit)[2]
> sig2.hat = sum(residuals(ar1.fit)^2)/ar1.fit$df.resid
> lrv.ar1 = sig2.hat/(1-rho.hat)^2
> as.numeric(lrv.ar1)
[1] 13.75 d Here lrvAR (yt ) = 13.75, and an estimate for SE(¯) is SE AR (¯) = 0.371. y y
The S+FinMetrics function asymp.var may be used to compute the nonparameteric Newey-West estimate lrvNW (yt ). The arguments expected by asymp.var are

88

3. Time Series Concepts

> args(asymp.var) function(x, bandwidth, window = "bartlett", na.rm = F) where x is a “timeSeries”, bandwidth sets the truncation lag MT in
(3.20) and window specifies the weight function. Newey and West suggest setting the bandwidth using the sample size dependent rule
MT = 4(T /100)2/9 which is equal to 4 in the present case. The Newey-West long-run variance estimate is then
> lrv.nw = asymp.var(y.ar1, bandwidth=4)
> lrv.nw
[1] 7.238 d and the Newey-West estimate of SE(¯) is SE N W (¯) = 0.269. y y

3.2.9 Variance Ratios

There has been considerable interest in testing the so-called random walk
(RW) model for log stock prices (see chapter 2 in Campbell, Lo and MacKinlay (1997) for an extensive review). The RW model for log prices pt has the form pt = µ + pt−1 + εt , t = 1, . . . , T where εt is a random error term. Using rt = ∆pt , the RW model may be rewritten as rt = µ + εt
Campbell, Lo and MacKinlay distinguish three forms of the random walk model: RW1 εt ∼ iid(0, σ 2 )
RW2 εt is an independent process (allows for heteroskedasticity)
RW3 εt is an uncorrelated process (allows for dependence in higher order moments) For asset returns, RW1 and RW2 are not very realistic and, therefore, most attention has been placed on testing the model RW3.
Some commonly used tests for RW3 are based on constructing variance ratios. To illustrate, consider the simple two-period variance ratio
V R(2) =

var(rt (2))
2 · var(rt )

3.2 Univariate Time Series

89

The numerator of the variance ratio is the variance of the two-period return, rt (2) = rt−1 + rt , and the deminator is two times the variance of the oneperiod return, rt . Under RW1, is easy to see that V R (2) = 1. If {rt } is an ergodic-stationary process then
V R(2) =
=

var(rt−1 ) + var(rt ) + 2 · cov(rt , rt−1 )
2 · var(rt )
2γ 0 + 2γ 1
= 1 + ρ1
2γ 0

There are three cases of interest depending on the value of ρ1 . If ρ1 = 0 then V R(2) = 1; if ρ1 > 1 then V R(2) > 1; if ρ1 < 1 then V R(2) < 1.
The general q−period variance ratio is
V R(q) =

var(rt (q)) q · var(rt )

(3.21)

where rt (q) = rt−q+1 + · · · + rt . Under RW1, V R(q) = 1. For ergodic stationary returns, some algebra shows that
¶
q
Xµ
k
V R(q) = 1 + 2 ·
1−
ρk q k=1

When the variance ratio is greater than one, returns are called mean averting due to the dominating presence of positive autocorrelations. When the variance ratio is less than one, returns are called mean reverting due to the dominating presence of negative autocorrelations. Using the Wold representation (3.6), it can be shown that lim V R(q) =

q→∞

σ 2 ψ(1)2 lrv(rt )
=
γ0 var(rt )

That is, as q becomes large the variance ratio approaches the ratio of the long-run variance to the short-run variance. Furthermore, Under RW2 and
RW3 it can be shown that V R(q) → 1 as q → ∞ provided
T
1X var(rt ) → σ 2 > 0
¯
T t=1

Test Statistics
Let {p0 , p1 , . . . , pT q } denote a sample of T q + 1 log prices, which produces a sample of T q one-period returns {r1 , . . . , rT q }. Lo and MacKinlay (1988,
1989) develop a number of test statistics for testing the random walk hypothesis based on the estimated variance ratio v ar(rt (q)) d d
V R(q) = q · v ar(rt ) d (3.22)

90

3. Time Series Concepts

The form of the statistic depends on the particular random walk model
(RW1, RW2 or RW3) assumed under the null hypothesis.
Under RW1, (3.22) is computed using σ 2 (q)
ˆ
d
V R(q) = σ2 ˆ

where σ2 ˆ

=

Tq
1 X
(rk − µ)2
ˆ
Tq k=1 σ 2 (q) =
ˆ

Tq
1 X
(rk (q) − qˆ )2 µ T q2 k=q µ =
ˆ

1
Tq

Tq
X

rk =

k=1

1
(pT q − p0 )
Tq

Lo and MacKinlay show that, under RW1, p A d T q(V R(q) − 1) ∼ N (0, 2(q − 1))

Therefore, the variance ratio test statistic
ˆ
ψ(q) =

µ

Tq
2(q − 1)

¶1/2

d
(V R(q) − 1)

(3.23)

has a limiting standard normal distribution under RW1.
Lo and MacKinlay also derive a modified version of (3.23) based on the following bias corrected estimates of σ 2 and σ 2 (q) :
Tq

σ2
¯

=

1 X
(rk − µ)2
ˆ
Tq − 1 k=1 σ 2 (q) =
¯

1 m Tq
X

k=q

(rk (q) − qˆ )2 µ µ
¶
q m = q(T q − q + 1) 1 −
Tq
¯ σ Defining V R(q) = σ 2 (q)/¯ 2 , the biased corrected version of (3.23) has the form µ
¶1/2
3T q 2
¯
ψ(q) =
(V R(q) − 1)
(3.24)
2(2q − 1)(q − 1) which has a limiting standard normal distribution under RW1.

3.2 Univariate Time Series

91

The variance ratio statistics (3.23) and (3.24) are not valid under the empirically relevant RW3 and RW3 models. For this model, Lo and MacKinlay derived the heteroskedasticity robust variance ratio statistic
ˆ
ψ ∗ (q) = Ω(q)−1/2 (V R(q) − 1)

(3.25)

where
ˆ
Ω(q) =

q−1
X µ 2(q − j) ¶

j

j=1

ˆj δ αjt
ˆ

=

ˆj δ PT q

α0t αjt
ˆ ˆ
´2
α0t
ˆ
j=1

t=j+1

³P
Tq

= (rt−j − rt−j−1 − µ)
ˆ

Under RW2 or RW3, Lo and MacKinlay show that (3.25) has a limiting stardard normal distribution.
Example 10 Testing the random walk hypothesis using variance ratios
The variance ratio statistics (3.23), (3.24) and (3.25) may be computed using the S+FinMetrics function varRatioTest. The arguments for varRatioTest are
> args(varRatioTest) function(x, n.periods, unbiased = T, hetero = F) where x is the log return series (which may contain more than one series) and n.periods denotes the number of periods q in the variance ratio. If unbiased=T and hetero=F the bias corrected test statistic (3.24) is computed. If unbiased=T and hetero=T then the heteroskedasticity robust statistic (3.25) is computed. The function varRatioTest returns an object of class “varRatioTest” for which there are print and plot methods.
Consider testing the model RW3 for the daily log closing prices of the
Dow Jones Industrial Average over the period 1/1/1960 through 1/1/1990.
To compute the variance ratio (3.21) and the heteroskedasticity robust test
(3.25) for q = 1, . . . , 60 use
> VR.djia = varRatioTest(djia[timeEvent("1/1/1960","1/1/1990"),
+
"close"], n.periods=60, unbiased=T, hetero=T)
> class(VR.djia)
[1] "varRatioTest"
> names(VR.djia)
[1] "varRatio" "std.err" "stat"
"hetero"
> VR.djia
Variance Ratio Test

92

3. Time Series Concepts

Variance Ratio Profile

1.2
1.0
0.4

0.6

0.8

Variance Ratio

1.4

1.6

1.8

close

0

10

20

30

40

50

60

Periods

FIGURE 3.12. Variance ratios for the daily log prices of the Dow Jones Industrial
Average.

Null Hypothesis: random walk with heteroskedastic errors
Variable: close var.ratio std.err
2
1.0403 0.06728
3
1.0183 0.10527
...
60
1.0312 0.36227

stat
0.5994871
0.1738146
0.0861747

* : significant at 5% level
** : significant at 1% level
None of the variance ratios are statistically different from unity at the 5% level. Figure 3.12 shows the results of the variance ratio tests based on plot method > plot(VR.djia)
The variance ratios computed for different values of q hover around unity, and the ± 2 × standard error bands indicate that the model RW3 is not rejected at the 5% level.

3.3 Univariate Nonstationary Time Series

93

Variance Ratio Profile

4

5

1.0

4

5

0.8
0.7
0.6
1.1

MSFT

0.9
0.8
1.10
0.90
0.80

HON

1.1

1.1

1.2

HD

0.8

0.9

0.9

1.0

1.0

0.90

JPM

1.00

1.1

GM

0.80

1.10

CAT

DD

1.00

1.00

C

0.80

0.80

0.8

0.90

0.90

0.70
1.0
0.9

0.95
0.85

3

XOM

1.0

1.00
0.80

JNJ

0.8

0.80
1.00

1.0
0.9
1.05 0.7

BA

2

0.9

1.1
0.9
0.8
1.10

MRK

0.9

0.90

0.8

GE

0.8
3

1

WMT

1.0

1.0

1.00

IP

0.75

0.8
0.7
2

5

0.90

0.90
0.80

0.80

INTC

0.7

0.7
1.15
1.05
0.95
1.0
0.9

1.0
0.9
0.8
1

4

1.0

1.1
0.9
0.8

1.00

1.00

MO

0.9

0.9
0.8

0.7
0.95
0.85
1.1
0.75

AXP

3

1.0

1.0
1.10

MMM

EK

AA

2

UTX

0.90

0.90
0.80
1.0

1.0
0.9
0.8

1

T

0.8

0.8
1.00

1.10

MCD

IBM

DIS

5

0.9

0.9

1.0
0.9
0.8
1.15
1.05
0.95
0.85

HWP

4

1.1

SBC

KO

1.05

Variance Ratio

3

1.0

1.1

PG

2

1.1

1

1

2

3

4

5

1

2

3

4

5

Periods

FIGURE 3.13. Variance ratio statistics for daily log prices on individual Dow
Jones index stocks.

The RW3 model appears to hold for the Dow Jones index. To test the
RW3 model for the top thirty stocks in the index individually, based on q = 1, . . . , 5, use
> VR.DJ30 = varRatioTest(DowJones30, n.periods=5, unbiased=T,
+
hetero=T)
> plot(VR.DJ30)
The results, illustrated in Figure 3.13, indicate that the RW3 model may not hold for some individual stocks.

3.3 Univariate Nonstationary Time Series
A univariate time series process {yt } is called nonstationary if it is not stationary. Since a stationary process has time invariant moments, a nonstationary process must have some time dependent moments. The most common forms of nonstationarity are caused by time dependence in the mean and variance.

94

3. Time Series Concepts

Trend Stationary Process
{yt } is a trend stationary process if it has the form yt = T Dt + xt where T Dt are deterministic trend terms (constant, trend, seasonal dummies etc) that depend on t and {xt } is stationary. The series yt is nonstationary because E[T Dt ] = T Dt which depends on t. Since xt is stationary, yt never deviates too far away from the deterministic trend T Dt . Hence, yt exhibits trend reversion. If T Dt were known, yt may be transformed to a stationary process by subtracting off the deterministic trend terms: xt = yt − T Dt
Example 11 Trend stationary AR(1)
A trend stationary AR(1) process with T Dt = µ + δt may be expressed in three equivalent ways yt = µ + δt + ut , ut = φut−1 + εt yt − µ − δt = φ(yt−1 − µ − δ(t − 1)) + εt yt = c + βt + φyt−1 + εt where |φ| < 1, c = µ(1 − φ) + δ, β = δ(1 − φ)t and εt ∼ W N (0, σ 2 ). Figure
3.14 shows T = 100 observations from a trend stationary AR(1) with µ = 1, δ = 0.25, φ = 0.75 and σ 2 = 1 created with the S-PLUS commands
>
>
+
>
>

set.seed(101)
y.tsar1 = 1 + 0.25*seq(100) + arima.sim(model=list(ar=0.75),n=100) tsplot(y.tsar1,ylab="y") abline(a=1,b=0.25) The simulated data show clear trend reversion.
Integrated Processes
{yt } is an integrated process of order 1, denoted yt ∼ I(1), if it has the form yt = yt−1 + ut

(3.26)

where ut is a stationary time series. Clearly, the first difference of yt is stationary ∆yt = ut
Because of the above property, I(1) processes are sometimes called difference stationary processes. Starting at y0 , by recursive substitution yt has

95

5

10

y

15

20

3.3 Univariate Nonstationary Time Series

0

20

40

60

80

100

FIGURE 3.14. Simulated trend stationary process.

the representation of an integrated sum of stationary innovations yt = y0 +

t
X

uj .

(3.27)

j=1

The integrated sum
T St . Notice that

Pt

j=1

uj is called a stochastic trend and is denoted
T St = T St−1 + ut

where T S0 = 0. In contrast to a deterministic trend, changes in a stochastic trend are not perfectly predictable.
Since the stationary process ut does not need to be differenced, it is called an integrated process of order zero and is denoted ut ∼ I(0). Recall, from the Wold representation (3.6) a stationary process has an infinite order moving average representation where the moving average weights decline to zero at a geometric rate. From (3.27) it is seen that an I(1) process has an infinite order moving average representation where all of the weights on the innovations are equal to 1.
If ut ∼ IW N (0, σ 2 ) in (3.26) then yt is called a random walk. In general, an I(1) process can have serially correlated and heteroskedastic innovations ut . If yt is a random walk and assuming y0 is fixed then it can be shown

96

3. Time Series Concepts

that γ0 γj ρj = σ2 t
= (t − j)σ 2 r t−j
=
t

which clearly shows that yt is nonstationary. Also, if t is large relative to j then ρj ≈ 1. Hence, for an I(1) process, the ACF does not decay at a geometric rate but at a linear rate as j increases.
An I(1) process with drift has the form yt = µ + yt−1 + ut , where ut ∼ I(0)
Starting at t = 0 an I(1) process with drift µ may be expressed as yt = y0 + µt +

t
X

ut

j=1

= T Dt + T St

so that it may be thought of as being composed of a deterministic linear
Pt
trend T Dt = y0 + µt as well as a stochastic trend T St = j=1 uj .
An I(d) process {yt } is one in which ∆d yt ∼ I(0). In finance and economics data series are rarely modeled as I(d) process with d > 2. Just as an I(1) process with drift contains a linear deterministic trend, an I(2) process with drift will contain a quadratic trend.
Example 12 Simulated I(1) processes
Consider the simulation of T = 100 observations from various I(1) processes where the innovations ut follow an AR(1) process ut = 0.75ut−1 + εt with εt ∼ GW N (0, 1).
>
>
>
>
>
>
+
+

set.seed(101)
u.ar1 = arima.sim(model=list(ar=0.75), n=100) y1 = cumsum(u.ar1) y1.d = 1 + 0.25*seq(100)+ y1 y2 = rep(0,100) for (i in 3:100) { y2[i] = 2*y2[i-1] - y2[i-2] + u.ar1[i]
}

The simulated data are illustrated in Figure 3.15 .
Example 13 Financial time series

3.4 Long Memory Time Series
I(1) process

-2

0

0

10

2

20

30

4

I(0) innovations

97

0

20

40

60

80

100

0

20

60

80

100

80

100

I(2) process

0

0

10

500

20

30

1000

40

50

1500

I(1) process with drift

40

0

20

40

60

80

100

0

20

40

60

FIGURE 3.15. Simulated I(d) processes for d = 0, 1 and 2.

Many financial time series are well characterized by I(1) processes. The leading example of an I(1) process with drift is the logarithm of an asset price. Common examples of I(1) processes without drifts are the logarithms of exchange rates, nominal interest rates, and inflation rates. Notice that if inflation is constructed as the the difference in the logarithm of a price index and is an I(1) process, then the logarithm of the price index is an
I(2) process. Examples of these data are illustrated in Figure 3.16. The exchange rate is the monthly log of the US/CA spot exchange rate taken from the S+FinMetrics “timeSeries” lexrates.dat, the asset price of the monthly S&P 500 index taken from the S+FinMetrics “timeSeries” object singleIndex.dat, the nominal interest rate is the 30 day T-bill rate taken from the S+FinMetrics “timeSeries” object rf.30day, and the monthly consumer price index is taken from the S+FinMetrics “timeSeries” object CPI.dat.

3.4 Long Memory Time Series
If a time series yt is I(0) then its ACF declines at a geometric rate. As a result, I(0) process have short memory since observations far apart in time are essentially independent. Conversely, if yt is I(1) then its ACF declines at a linear rate and observations far apart in time are not independent. In

98

3. Time Series Concepts

Log S&P 500 index

5.8

-0.35

6.2

6.6

-0.15

7.0

0.00

Log US/CA spot exchange rate

1976

1980

1984

1988

1992

1996

1990

1992

1996

1998

2000

Log of US CPI

2.5

0.000

3.0

3.5

4.0

0.008

4.5

5.0

Nominal 30 day T-bill rate

1994

1930

1950

1970

1990

1915

1935

1955

1975

1995

FIGURE 3.16. Monthly financial time series.

between I(0) and I(1) processes are so-called fractionally integrated I(d) process where 0 < d < 1. The ACF for a fractionally integrated processes declines at a polynomial (hyperbolic) rate, which implies that observations far apart in time may exhibit weak but non-zero correlation. This weak correlation between observations far apart is often referred to as long memory.
A fractionally integrated white noise process yt has the form
(1 − L)d yt = εt , εt ∼ W N (0, σ 2 )

(3.28)

d

where (1 − L) has the binomial series expansion representation (valid for any d > −1)
∞
Xµ d ¶ d (−L)k
(1 − L) = k k=0

d(d − 1) 2 d(d − 1)(d − 2) 3
L −
L + ···
2!
3!
If d = 1 then yt is a random walk and if d = 0 then yt is white noise. For
0 < d < 1 it can be shown that
= 1 − dL +

ρk ∝ k 2d−1 as k → ∞ so that the ACF for yt declines hyperbolically to zero at a speed that depends on d. Further, it can be shown yt is stationary and ergodic for 0 < d < 0.5 and that the variance of yt is infinite for 0.5 ≤ d < 1.

99

-2

0

2

4

3.4 Long Memory Time Series

0

100

200

300

400

500

40

50

ACF

0.0

0.2

0.4

0.6

0.8

1.0

Series : y.fwn

0

10

20

30
Lag

FIGURE 3.17. Simulated values from a fractional white noise process with d = 0.3 and σ = 1.

Example 14 Simulated fractional white noise
The S+FinMetrics function simulate.FARIMA may be used to generate simulated values from a fractional white noise process. To simulate 500 observations from (3.28) with d = 0.3 and σ 2 = 1 use
> set.seed(394)
> y.fwn = simulate.FARIMA(list(d=0.3), 500)
Figure 3.17 shows the simulated data along with the sample ACF created using > par(mfrow=c(2,1))
> tsplot(y.fwn)
> tmp = acf(y.fwn,lag.max=50)
Notice how the sample ACF slowly decays to zero.
A fractionally integrated process with stationary and ergodic ARMA(p, q) errors (1 − L)d yt = ut , ut ∼ ARM A(p, q) is called an autoregressive fractionally integrated moving average (ARFIMA) process. The modeling of long memory process is described in detail in
Chapter 8.

100

3. Time Series Concepts

0.6
0.4
0.0

0.2

ACF

0.8

1.0

Series : msft.aret

0

20

40

60

80

100

40

50

Lag

0.4
-0.2

0.2

ACF

0.6

0.8

1.0

Series : uscn.id

0

10

20

30
Lag

FIGURE 3.18. SACFs for the absolute value of daily returns on Microsoft and the monthly 30-day interest rate differential between U.S. bonds and Canadian bonds. Example 15 Long memory in financial time series
Long memory behavior has been observed in certain types of financial time series. Ding, Granger and Engle (1993) find evidence of long memory in the absolute value of daily stock returns. Baillie and Bollerslev (1994) find evidence for long memory in the monthly interest rate differentials between short term U.S. government bonds and short term foreign government bonds. To illustrate, consider the absolute values of the daily returns on Microsoft over the 10 year period 1/2/1991 - 1/2/2001 taken from the
S+FinMetrics “timeSeries” DowJones30
> msft.aret = abs(getReturns(DowJones30[,"MSFT"]))
Consider also the monthly US/CA 30-day interest rate differential over the period February 1976 through June 1996 in the “timeSeries” uscn.id constructed earlier and taken from the S+FinMetrics “timeSeries” object lexrates.dat. Figure 3.18 shows the SACFs these series create by
> par(mfrow=c(2,1))
> tmp = acf(msft.aret, lag.max=100)
> tmp = acf(uscn.id, lag.max=50)
For the absolute return series, notice the large number of small but apparently significant autocorrelations at very long lags. This is indicative of

3.5 Multivariate Time Series

101

long memory. For the interest rate differential series, the ACF appears to decay fairly quickly, so the evidence for long memory is not as strong.

3.5 Multivariate Time Series
Consider n time series variables {y1t }, . . . , {ynt }. A multivariate time series is the (n × 1) vector time series {Yt } where the ith row of {Yt } is
{yit }. That is, for any time t, Yt = (y1t , . . . , ynt )0 . Multivariate time series analysis is used when one wants to model and explain the interactions and co-movements among a group of time series variables. In finance, multivariate time series analysis is used to model systems of asset returns, asset prices and exchange rates, the term structure of interest rates, asset returns/prices, and economic variables etc. Many of the time series concepts described previously for univariate time series carry over to multivariate time series in a natural way. Additionally, there are some important time series concepts that are particular to multivariate time series. The following sections give the details of these extensions and provide examples using
S-PLUS and S+FinMetrics.

3.5.1 Stationary and Ergodic Multivariate Time Series
A multivariate time series Yt is covariance stationary and ergodic if all of its component time series are stationary and ergodic. The mean of Yt is defined as the (n × 1) vector
E[Yt ] = µ = (µ1 , . . . , µn )0 where µi = E[yit ] for i = 1, . . . , n. The variance/covariance matrix of Yt is the (n × n) matrix var(Yt ) = Γ0 = E[(Yt −µ)(Yt −µ)0 ]

var(y1t ) cov(y1t , y2t )
 cov(y2t , y1t ) var(y2t )

= 
.
.
.
.

.
.

···
···
..
.

cov(ynt , y1t ) cov(ynt , y2t ) · · ·

cov(y1t , ynt ) cov(y2t , ynt )
.
.
.
var(ynt )







The matrix Γ0 has elements γ 0 = cov(yit , yjt ). The correlation matrix of ij Yt is the (n × n) matrix corr(Yt ) = R0 = D−1 Γ0 D−1 where D is an (n×n) diagonal matrix with j th diagonal element (γ 0 )1/2 = jj SD(yjt ). The parameters µ, Γ0 and R0 are estimated from data (Y1 , . . . , YT )

102

3. Time Series Concepts

using the sample moments
T
1X
Yt
T t=1

¯
Y

=

ˆ
Γ0

=

ˆ
R0

ˆ ˆ ˆ
= D−1 Γ0 D−1

T
1X
¯
¯
(Yt −Y)(Yt −Y)0
T t=1

where D is the (n×n) diagonal matrix with the sample standard deviations
ˆ
of yjt along the diagonal. In order for the sample variance matrix Γ0 and
ˆ 0 to be positive definite, the sample size T must be correlation matrix R greater than the number of component time series n.
Example 16 System of asset returns
The S+FinMetrics “timeSeries” object DowJones30 contains daily closing prices on the 30 assets in the Dow Jones index. An example of a stationary and ergodic multivariate time series is the continuously compounded returns on the first four assets in this index:
> Y = getReturns(DowJones30[,1:4],type="continuous")
> colIds(Y)
[1] "AA" "AXP" "T"
"BA"
The S-PLUS function colMeans may be used to efficiently compute the mean vector of Y
> colMeans(seriesData(Y))
AA
AXP
T
BA
0.0006661 0.0009478 -0.00002873 0.0004108
The function colMeans does not have a method for “timeSeries” objects so the extractor function seriesData is used to extract the data slot of the variable Y. The S-PLUS functions var and cor, which do have methods
ˆ
ˆ for “timeSeries” objects, may be used to compute Γ0 and R0
> var(Y)
AA
AXP
T
AA 0.00041096 0.00009260 0.00005040
AXP 0.00009260 0.00044336 0.00008947
T 0.00005040 0.00008947 0.00040441
BA 0.00007301 0.00009546 0.00004548
> cor(Y)
AA
AXP
T
BA
AA 1.0000 0.2169 0.1236 0.1877
AXP 0.2169 1.0000 0.2113 0.2362
T 0.1236 0.2113 1.0000 0.1179

BA
0.00007301
0.00009546
0.00004548
0.00036829

3.5 Multivariate Time Series

103

BA 0.1877 0.2362 0.1179 1.0000
If only the variances or standard deviations of Yt are needed the S-PLUS functions colVars and colStdevs may be used
> colVars(seriesData(Y))
AA
AXP
T
BA
0.000411 0.0004434 0.0004044 0.0003683
> colStdevs(seriesData(Y))
AA
AXP
T
BA
0.020272 0.021056 0.02011 0.019191

Cross Covariance and Correlation Matrices
For a univariate time series yt the autocovariances γ k and autocorrelations ρk summarize the linear time dependence in the data. With a multivariate time series Yt each component has autocovariances and autocorrelations but there are also cross lead-lag covariances and correlations between all possible pairs of components. The autocovariances and autocorrelations of yjt for j = 1, . . . , n are defined as γk jj

= cov(yjt , yjt−k ),

ρk jj = corr(yjt , yjt−k ) =

γk jj γ0 jj and these are symmetric in k: γ k = γ −k , ρk = ρ−k . The cross lag covarijj jj jj jj ances and cross lag correlations between yit and yjt are defined as γk ij ρk ij

= cov(yit , yjt−k ), γk ij
= corr(yjt , yjt−k ) = q
0 γ0 γ ii jj

and they are not necessarily symmetric in k. In general, γ k = cov(yit , yjt−k ) 6= cov(yit , yjt+k ) = cov(yjt , yit−k ) = γ −k ij ij
If γ k 6= 0 for some k > 0 then yjt is said to lead yit . Similarly, if γ −k 6= 0 ij ij for some k > 0 then yit is said to lead yjt . It is possible that yit leads yjt and vice-versa. In this case, there is said to be feedback between the two series. 104

3. Time Series Concepts

All of the lag k cross covariances and correlations are summarized in the
(n × n) lag k cross covariance and lag k cross correlation matrices
Γk

Rk

= E[(Yt −µ)(Yt−k −µ)0 ]

cov(y1t , y1t−k ) cov(y1t , y2t−k )
 cov(y2t , y1t−k ) cov(y2t , y2t−k )

= 
.
.
.
.

.
.

···
···
..
.

cov(y1t , ynt−k ) cov(y2t , ynt−k )
.
.
.

cov(ynt , y1t−k ) cov(ynt , y2t−k ) · · ·

cov(ynt , ynt−k )

−1

= D

−1

Γk D







The matrices Γk and Rk are not symmetric in k but it is easy to show that
Γ−k = Γ0 and R−k = R0 . The matrices Γk and Rk are estimated from data k k
(Y1 , . . . , YT ) using
T
1 X
¯
¯
(Yt −Y)(Yt−k −Y)0
T

ˆ
Γk

=

ˆ
Rk

ˆ ˆ ˆ
= D−1 Γk D−1

t=k+1

Example 17 Lead-lag covariances and correlations among asset returns
Consider computing the cross lag covariances and correlations for k =
0, . . . , 5 between the first two Dow Jones 30 asset returns in the “timeSeries”
Y. These covariances and correlations may be computed using the S-PLUS function acf
> Ghat = acf(Y[,1:2],lag.max=5,type="covariance",plot=F)
> Rhat = acf(Y[,1:2],lag.max=5,plot=F)
Ghat and Rhat are objects of class “acf” for which there is only a print method. For example, the estimated cross lag autocorrelations are
> Rhat
Call: acf(x = Y[, 1:2], lag.max = 5, plot = F)
Autocorrelation matrix: lag AA.AA AA.AXP AXP.AXP
1
0 1.0000 0.2169 1.0000
2
1 0.0182 0.0604 -0.0101
3
2 -0.0556 -0.0080 -0.0710
4
3 0.0145 -0.0203 -0.0152
5
4 -0.0639 0.0090 -0.0235
6
5 0.0142 -0.0056 -0.0169

1

lag
0

AXP.AA
0.2169

3.5 Multivariate Time Series

105

Multivariate Series : Y[, 1:2]
AA and AXP

0.0

0.0

0.2

0.05

ACF
0.4

0.10

0.6

0.15

0.8

0.20

1.0

AA

0

1

2

3

4

5

0

1

2

4

5

3

4

5

AXP

0.0

0.0

0.2

0.05

0.4

ACF
0.10

0.6

0.15

0.8

0.20

1.0

AXP and AA

3

-5

-4

-3

Lag

-2

-1

0

0

1

2

Lag

FIGURE 3.19. Cross lag correlations between the first two Dow Jones 30 asset returns. 2
3
4
5
6

-1
-2
-3
-4
-5

-0.0015
-0.0187
-0.0087
-0.0233
0.0003

The function acf.plot may be used to plot the cross lag covariances and correlations produced by acf.
> acf.plot(Rhat)
ˆ
ˆ
Figure 3.19 shows these cross lag correlations. The matrices Γk and Rk may be extracted from acf component of Ghat and Rhat, respectively. For example, > Ghat$acf[1,,]
[,1]
[,2]
[1,] 0.00041079 0.00009256
[2,] 0.00009256 0.00044318
> Rhat$acf[1,,]
[,1]
[,2]
[1,] 1.0000 0.2169
[2,] 0.2169 1.0000
> Ghat$acf[2,,]

106

3. Time Series Concepts

[,1]
[,2]
[1,] 7.488e-006 2.578e-005
[2,] -6.537e-007 -4.486e-006
> Rhat$acf[2,,]
[,1]
[,2]
[1,] 0.018229 0.06043
[2,] -0.001532 -0.01012
ˆ ˆ ˆ
ˆ
extracts Γ1 , R1 , Γ2 and R2 .

3.5.2 Multivariate Wold Representation
Any (n × 1) covariance stationary multivariate time series Yt has a Wold or linear process representation of the form
Yt

= µ + εt +Ψ1 εt−1 +Ψ2 εt−2 + · · ·
∞
X
= µ+
Ψk εt−k

(3.29)

k=0

where Ψ0 = In and εt is a multivariate white noise process with mean zero and variance matrix E[εt ε0 ] = Σ. In (3.29), Ψk is an (n × n) matrix with t (i, j)th element ψ k . In lag operator notation, the Wold form is ij = µ + Ψ(L)εt
∞
X
Ψ(L) =
Ψk Lk
Yt

k=0

The moments of Yt are given by
E[Yt ] = µ
∞
X var(Yt ) =
Ψk ΣΨ0 k k=0

VAR Models
The most popular multivariate time series model is the vector autoregressive
(VAR) model. The VAR model is a multivariate extension of the univariate autoregressive model. For example, a bivariate VAR(1) model has the form µ ¶ µ
¶ µ 1
¶µ
¶ µ
¶
y1t c1 π 11 π 1 y1t−1 ε1t
12
=
+
+ y2t c2 π1 π1 y2t−1 ε2t
21
22 or y1t y2t = c1 + π 1 y1t−1 + π1 y2t−1 + ε1t
11
12
= c2 + π 1 y1t−1 + π1 y2t−1 + ε2t
21
22

3.5 Multivariate Time Series

where

µ

ε1t ε2t ¶

∼ iid

µµ

0
0

¶ µ σ 11
,
σ 12

σ 12 σ 22

107

¶¶

In the equations for y1 and y2 , the lagged values of both y1 and y2 are present. The general VAR(p) model for Yt = (y1t , y2t , . . . , ynt )0 has the form
Yt = c + Π1 Yt−1 +Π2 Yt−2 + · · · + Πp Yt−p + εt , t = 1, . . . , T

(3.30)

where Πi are (n × n) coefficient matrices and εt is an (n × 1) unobservable zero mean white noise vector process with covariance matrix Σ. VAR models are capable of capturing much of the complicated dynamics observed in stationary multivariate time series. Details about estimation, inference, and forecasting with VAR models are given in chapter eleven.

3.5.3 Long Run Variance
Let Yt be an (n × 1) stationary and ergodic multivariate time series with
E[Yt ] = µ. Anderson’s central limit theorem for stationary and ergodic process states


∞
X
√
d
¯
T (Y − µ) → N 0,
Γj  j=−∞ or




∞
X
1
¯ A
Y ∼ N µ,
Γj 
T j=−∞

Hence, the long-run variance of Yt is T times the asymptotic variance of
¯
Y:
∞
X
¯
lrv(Yt ) = T · avar(Y) =
Γj
j=−∞

Since Γ−j = Γ0 , lrv(Yt ) may be alternatively expressed as j lrv(Yt ) = Γ0 +

∞
X
(Γj +Γ0 ) j j=1

Using the Wold representation of Yt it can be shown that lrv(Yt ) = Ψ(1)ΣΨ(1)0 where Ψ(1) =

P∞

k=0

Ψk .

108

3. Time Series Concepts

VAR Estimate of the Long-Run Variance
The Wold representation (3.29) may be approximated by high order VAR(p∗ ) model Yt = c + Φ1 Yt−1 + · · · + Φp∗ Yt−p∗ +εt where the lag length p∗ is chosen such p∗ = O(T 1/3 ). This gives rise to the autoregressive long-run variance matrix estimate c ˆ
ˆˆ 0 lrv AR (Yt ) = Ψ(1)ΣΨ(1)
ˆ
ˆ
ˆ
Ψ(1) = (In − Φ1 − · · · − Φp )−1
ˆ
Σ =

T
1X
ˆtˆ0 ε εt
T t=1

(3.31)
(3.32)
(3.33)

ˆ where Φk (k = 1, . . . , p∗ ) are estimates of the VAR parameter matrices.
Non-parametric Estimate of the Long-Run Variance
A consistent estimate of lrv(Yt ) may be computed using non-parametric methods. A popular estimator is the Newey-West weighted autocovariance estimator MT
³
´
X
c
ˆ
ˆ
ˆj
lrv NW (Yt ) = Γ0 + wj,T · Γj + Γ0
(3.34)
j=1

where wj,T are weights which sum to unity and MT is a truncation lag parameter that satisfies MT = O(T 1/3 ).

Example 18 Newey-West estimate of long-run variance matrix for stock returns The S+FinMetrics function asymp.var may be used to compute the
Newey-West long-run variance estimate (3.34) for a multivariate time series.
The long-run variance matrix for the first four Dow Jones assets in the
“timeSeries” Y is
> M.T = floor(4*(nrow(Y)/100)^(2/9))
> lrv.nw = asymp.var(Y,bandwidth=M.T)
> lrv.nw
AA
AXP
T
BA
AA 0.00037313 0.00008526 3.754e-005 6.685e-005
AXP 0.00008526 0.00034957 7.937e-005 1.051e-004
T 0.00003754 0.00007937 3.707e-004 7.415e-006
BA 0.00006685 0.00010506 7.415e-006 3.087e-004

3.6 References

109

3.6 References
[1] Alexander, C. (2001). Market Models. A Guide to Financial Data
Analysis. John Wiley & Sons, Chichester, UK.
[2] Baille, R.T. and T. Bollerslev (1994). “The Long Memory of the Forward Premium,” Journal of International Money and Finance,
13, 555-571.
[3] Box, G.E.P. and G.M. Jenkins (1976). Time Series Analysis,
Forecasting and Control. Revised Edition. Holden Day, San Francisco.
[4] Campbell, J.Y., A.W. Lo, A.C. MacKinlay (1997). The Econometrics of Financial Markets. Princeton University Press, New Jersey.
[5] Box, G.E.P., and D.A. Pierce (1970). “Distribution of Residual
Autocorrelations in Autoregressive-integrated Moving Average Time
Series Models,” Journal of the American Statistical Association, 65,
1509-1526.
[6] Chan, N.H. (2002). Time Series: Applicatios to Finance. John Wiley
& Sons, New York.
[7] Ding, Z., C.W.J. Granger and R.F. Engle (1993). “A Long
Memory Property of Stock Returns and a New Model,” Journal of
Empirical Finance, 1, 83-106.
[8] Engle, R.F. (1982). “Autoregressive Conditional Heteroskedasticity with Estimates of the Variance of United Kingdom Inflations,”
Econometrica, 50, 987-1097.
[9] Fuller, W.A. (1996). Introduction to Statistical Time Series, Second Edition. John Wiley & Sons, New York.
[10] Gourieroux, C and J. Jasiak (2001). Financial Econometrics.
Princeton University Press, New Jersey.
[11] Granger, C.W.J. and M.J. Morris (1976). “Time Series Modeling and Interpretation,” Journal of the Royal Statistical Society,
Series A, 139, 246-257.
[12] Hamilton, J.D. (1994). Time Series Analysis. Princeton University
Press, New Jersey.
[13] Harvey, A.C. (1993). Time Series Models, Second Edition. MIT
Press, Massachusetts.
[14] Hayashi, F. (2000). Econometrics. Princeton University Press, New
Jersey.

110

3. Time Series Concepts

[15] Ljung, T. and G.E.P. Box (1979). “The Likelihood Function for a Stationary Autoregressive Moving Average Process,” Biometrika,
66, 265-270.
[16] Mills, T.C. (1999). The Econometric Modeling of Financial Time
Series, Second Edition. Cambridge University Press, Cambridge, UK.
[17] Newey, W.K. and K.D. West (1987). “A Simple Positive Semidefinite Heteroskedasticity and Autocorrelation Consistent Covariance
Matrix,” Econometrica, 55, 703-708.
[18] Tsay, R.S. (2001). Analysis of Financial Time Series, John Wiley
& Sons, New York.
[19] Venables, W.N. and B.D. Ripley (1999). Modern Applied Statistics with S-PLUS, Third Edition. Springer-Verlag, New York.
[20] White, H. (1984). Asymptotic Theory for Econometrians. Academic
Press, San Diego.