Time Series Analysis II

Univariate ARIMA Models

Table of Contents

The Theoretical Model

• MA(q)
• AR(p)
• ARMA(p,q)

The Empirical Model

• Model Identification
• Model Estimation
• Forecasting

Extensions

• Intervention Analysis
• Seasonal ARMA and Mixed Model
• Transfer Function or ARMAX Model
• Multivariate Vector Autocorrelation or VAR Model

Readings

• W. H. Green, Chapter 20.
• K.-P. Lin, Chapter 15.

The General Model

or

Y_t = m + (1+y₁B+y₂B²+...)e_t = m + y(B)e_t

where

y(B) = 1 + y₁B + y₂B² + ..., and
e_t ~ ii(0,s²), t = 1,2,...,N.

Stationarity requirement for the process:

Mean	m = E(Yt) < ¥
Variance	g0 = s2åi=0,...,¥yi2 < ¥
Autocovariance	gj = s2åi=0,...,¥yiyj+i < ¥
Autocorrelation	fj = gj / g0

The plot of aucorrelation coefficients at each lag j = 0,1,2,... is called autocorrelation function. In general, the autocorrelation function of a stationary time series falls to 0 quickly. The autocorrelation function of a non-stationary time series does not converge to zero. However, by differencing the time series, the differenced series may become stationary.

A linear stochastic process called an integrated process of order d, I(d), if the d-th difference of the series is stationary. Note that d is the lowest number of differences required for the resulting series to be stationary. That is, Y_t ~ I(d) if D^dY_t is stationary, where

DY_t = Y_t-Y_t-1
D²Y_t = DY_t-DY_t-1
...

Moving Average Process of Order q: MA(q)

yj =	-qj, j=1,2,...q
	0 for j > q

then

Yt	= m - q1et-1 - q2et-2 - ... - qqet-q + et
	= m + (1 - q1B - q2B2 - ... - qqBq)et
	= m + q(B) et

where

q(B) = 1 - q₁B - q₂B² - ... - q_qB^q, and
e_t ~ ii(0,s²), t = 1,2,...,N.

• Stability Condition
  The roots (solutions) of the q-th order polynomial function of B, q(B), must lie outside the unit circle.

• Mean
  E(Y_t) = m

• Variance
  g₀ = Var(Y_t) = E[(Y_t-E(Y_t))²]
  Let y_t = Y_t-E(Y_t) = Y_t-m = q(B)e_t. Then,
  g₀ = E(y_t²) = E[(q(B)e_t)²] = (1+q₁²+q²+...+q_q²)s²

• Autocovariance
  g₁ = E(y_ty_t-1) = (-q₁+q₁q₂+q₂q₃+...+q_q-1q_q)s²
  g₂ = E(y_ty_t-2) = (-q₂+q₁q₃+q₂q₄+...+q_q-2q_q)s²
  ...
  g_q = E(y_ty_t-q) = -q_qs²

  or

  gj =	(-qj+q1qj+1+q2qj+2+...+qq-jqq)s2, j=1,2,...q
  	0 otherwise

• Autocorrelation (Normalized Autocovariance)
  f_j = g_j / g₀, j = 1,2,...

  or

  fj =	-qj+q1qj+1+q2qj+2+...+qq-jqq 1+q12+q22+...+qq2	, j = 1,2,...q
  	0	otherwise

  Plot of autocorrelation coefficients of a stationary MA(q) process indicates a finite memory pattern up to the q-th lag. Beyond that the autocorrelation function is zero-valued.

Examples

• MA(1): Y_t = m - q₁e_t-1 + e_t
• MA(2): Y_t = m - q₁e_t-1 - q₂e_t-2 + e_t

Autoregressive Process of Order p: AR(p)

Y_t = d + r₁Y_t-1 + r₂Y_t-2 + r₃Y_t-3 + ... e_t

Consider only the case with finite number of autoregressive lags. That is,

Y_t = d + r₁Y_t-1 + r₂Y_t-2 + ... + r_pY_t-p + e_t

or

r(B)Y_t = d + e_t,

where

r(B) = 1 - r₁B - r₂B² - ... - r_pB^p, and
e_t ~ ii(0,s²), t = 1,2,...,N.

• Stability Condition
  The roots (solutions) of the p-th order polynomial function of B, r(B), must lie outside the unit circle.

• Mean:
  m = E(Y_t) = d / (1-r₁-r₂-...-r_p)

• Variance:
  g₀ = Var(Y_t) = E[(Y_t-E(Y_t))²]
  Let y_t = Y_t-E(Y_t) = Y_t-m, then the AR(p) process can be written as:
  y_t = r₁y_t-1 + r₂y_t-2 + ... + r_py_t-p + e_t. Therefore

  g0 = E(yt2)	= E[(r1yt-1 + r2yt-2 + ... + rpyt-p + et)2]
  	= r1g1 + r2g2 + ... + rpgp + s2

• Autocovariance
  g₁ = E(y_ty_t-1) = r₁g₀ + r₂g₁ + r₃g₂ + ... + r_pg_p-1
  g₂ = E(y_ty_t-2) = r₁g₁ + r₂g₀ + r₃g₁ + ... + r_pg_p-2
  ...
  g_p = E(y_ty_t-p) = r₁g_p-1 + r₂g_p-2 + r₃g_p-3 + ... + r_pg₀
  ...
  g_j = E(y_ty_t-j) = r₁g_j-1 + r₂g_j-2 + r₃g_j-3 + ... + r_pg_j-p, for j > p

• Autocorrelation (Normalized Autocovaiance):
  f_j = g_j / g₀, j = 1,2,...

  That is the Yule-Walker Equations of AR(p) Process:
  f₁ = r₁ + r₂f₁ + r₃f₂ + ... + r_pf_p-1
  f₂ = r₁f₁ + r₂ + r₃f₁ + ... + r_pf_p-2
  ...
  f_p = r₁f_p-1 + r₂f_p-2 + r₃f_p-3 + ... + r_p
  ...
  f_j = r₁f_j-1 + r₂f_j-2 + r₃f_j-3 + ... + r_pf_j-p, for j > p

  Plot of autocorrelation coefficients f_j (j=1,2,...) of a stationary AR(p) process indicates an infinite but decay memory pattern.

Examples

• AR(1): Y_t = d + r₁Y_t-1 + e_t
  The Yule-Walker Equations:
  f₁ = r₁
  f_j = r₁f_j-1, for j > 1.
  That is,
  f_j = r₁^j, j=1,2,...

• AR(2): Y_t = d + r₁Y_t-1 + r₂Y_t-2 + e_t
  The Yule-Walker Equations:
  f₁ = r₁ + r₂f₁ (or f₁ = r₁/(1-r₂))
  f₂ = r₁f₁ + r₂
  ...
  f_j = r₁f_j-1 + r₂f_j-2, j > 2.

Mixed Autoregressive and Moving Average Process: ARMA(p,q)

or

r(B)Y_t = d + q(B)e_t,

where

r(B) = 1 - r₁B - r₂B² - ... - r_pB^p,
q(B) = 1 - q₁B - q₂B² - ... - q_qB^q, and
e_t ~ ii(0,s²), t = 1,2,...,N.

• Stability Condition
  The roots (solutions) of the p-th order polynomial function of B for the AR process, r(B), and the q-th order polynomial function of B for the MA process, q(B), must lie outside the unit circle, respectively.

• Mean:
  m = E(Y_t) = d / (1-r₁-r₂-...-r_p)

• Variance:
  g₀ = Var(Y_t) = E[(Y_t-E(Y_t))²]
  Let y_t = Y_t-E(Y_t) = Y_t-m, then the ARMA(p,q) process can be written as:
  y_t = r₁y_t-1 + r₂y_t-2 + ... + r_py_t-p - q₁y_t-1 - q₂y_t-2 - ... + q_qy_t-q + e_t.

  Therefore,

  g₀ = E(y_t²) = å_i=1,...,p r_ig_i + [1-å_i=1,...,q q_i(r_i-q_i)]s². That is,
  

  s2 =	g0-åi=1,...,p rigi 1-åi=1,...,q qi(ri-qi)

• Autocovariance
  g₁ = E(y_ty_t-1) = å_i=1,...,p r_ig_|i-1| - [q₁+å_i=1,...,q-1 q_i+1(r_i-q_i)]s².
  g₂ = E(y_ty_t-2) = å_i=1,...,p r_ig_|i-2| - [q₂+å_i=1,...,q-2 q_i+2(r_i-q_i)]s².
  ...

  or,

  for j = 1,2,...,q
  g_j = å_i=1,...,p r_ig_|i-j| - [q_j+å_i=1,...,q-j q_i+j(r_i-q_i)]s².
  
  and, for j > p
  g_j = å_i=1,...,p r_ig_|i-j|

• Autocorrelation (Normalized Autocovaiance):
  f_j = g_j / g₀, j = 1,2,...

  or,

  for j = 1,2,...,q

  fj =

  åi=1,...,p rif|i-j| -

  [qj+åi=1,...,q-j qi+j(ri-qi)] (1-åi=1,...,p rifi) 1-åi=1,...,q qi(ri-qi)

  and for, j > p
  f_j = å_i=1,...,p r_if_|i-j|

  Plot of autocorrelation coefficients f_j (j=1,2,...) of a stationary mixed ARMA(p,q) process reflects the combination of infinite memory AR and finite memory MA processes. However, after q lags, only the AR process continues.

Examples

• ARMA(1,1): Y_t = d + r₁Y_t-1 - q₁e_t-1 + e_t

  f1 =	(1-q1r1)(r1-q1) 1-2q1r1+q12

  f_j = r_jf_j-1, for j > 1.

• ARMA(2,2): Y_t = d + r₁Y_t-1 + r₂Y_t-2 + ... - q₁e_t-1 - q₂e_t-2 - ... + e_t

Partial Autocorrelation Function

f₁₁ = f₁ (this is r₁ of AR(1))
f₂₂ = (f₂-f₁²) / (1-f₁²) (this is r₂ of AR(2)),
and for additional lags j = 3,4,...:

fjj =

fj-åk=1,...,j-1 fj-1,kfj-k 1-åk=1,...,j-1 fj-1,kfk

where f_jk = f_j-1,k - f_jjf_j-1,j-k for k = 1,2,...,j-1.

Plot of partial autocorrelation can reveal the correct order of AR(p) process. In other words, For AR(p), f_jj = 0 for j > p. For MA(q), f_jj is non-zeros for all j and exhibits a geometrically decaying pattern. For ARMA(p,q), the decay of partial autocorrelations f_jj starts after the p-th lag.

The Empirical Model

• m = å_t=1,...,NY_t/N
• g₀ = å_t=1,...,N (Y_t-m)²/N
• g_j = å_t=1,...,N (Y_t-m)(Y_t+j-m)/N
• f_j = g_j/g₀, j = 1,2,...

The structural identification of a time series includes (1) the minimum order d of differencing required on the sample to achieve stationarity; (2) the appropriate order q of a moving-average process; and (3) the appropriate order p of an autoregressive process for the stationary time series. First of all, a rapidly decline pattern of sample autocorrelation plot or correlogram is needed to ensure a stationary time series for further identification and analysis.

Model Identification

• Autocorrelation Function (for identifying q)

  Bartlett Test

  Testing H₀: f_j = 0 for each j > q (no autocorrelation at each lag j longer than q) based on the Bartlett distribution of the estimated f_j.
  That is, the estimated f_j ~ normal(0,ÖVar(f_j)) approximately, where
  Var(f_j) = 1/N (1 + 2 å_j=1,...,qf_j²).

  Box-Pierece Test and Ljung-Box Test

  Testing H₀: f₁ =... = f_k =0 (zero autocorrelation coefficients up to some k lags) based on Box-Pierece Q or Ljung-Box Q' Statistic defined as follows:

  Q = N å_j=1,...,kf_j²
  Q' = N(N+2) å_j=1,...,kf_j²/(N-j)

  Both Q and Q' ~ Chi-Square(k).

• Partial Autocorrelation Function (for identifying p)

  Let f_jj be the partial autocorrelation coefficient at the j-th lag. That is, f_jj = r_j obtained from the autoregressive regression of the AR(j) model:

  Y_t = d + r₁Y_t-1 + r₂Y_t-2 + ... + r_jY_t-j + e_t

  If the sample series follows a AR(p) process, then the autoregressive coefficient for each lag j longer than p must be zero.

  Testing H₀: f_jj = 0 for each j > p based on the approximated distribution for the estimated f_jj ~ normal(0,1/ÖN). Alternatively, the standard error of f_jj can be obtained from the corresponding estimated autoregressive regression equation for each lag.

Model Estimation

• Pre-sample Data Initialization

  In order to use all N data observations, initialization may be needed for the following:

  • Y₀, Y_-1, ..., Y_-p+1 with E(Y_t) = d / (1-r₁-...-r_p)
    Note: å_j=1,...,p r_j ¹ 1.
  • e₀, e_-1, ..., e_-q+1 with E(e_t) = 0.

• Conditional Maximum Likelihood Estimation

  The model may be written in the "inverted" form as

  q(B)^-1(-d+r(B)Y_t) = e_t

  where e_t ~ ii(0,s²). Conditional to the historical information (Y_N, ..., Y₁), and data initialization (Y₀, ..., Y_-p+1), (e₀, ..., e_-q+1), the sum-of-squares is defined by

  S = å_t=1,2,...,Ne_t²

  Assuming e_t ~ nii(0,s²) for each observation i, the concentrated log-likelihood function is

  ll = -N/2 (1+ln(2p)-ln(N)+ln(S))

  The conditional maximum likelihood estimators of rs, qs, and d are obtained by maximizing the nonlinear function ll. A set of initial values for the parameters rs and qs are needed to start the iteration of nonlinear model estimation.

• Diagnosis Checking

  Further identification on the estimated residuals, with N-(p+q+1) of observations as the estimated model is an ARMA(p,q) process.

  • Residual analysis for normality and correlation
  • Autocorrelation correlation coefficients on the residuals: additional MA orders greater than q?
  • Partial autocorrelation coefficients on the residuals: additional AR orders greater than p?
  • Box-Pierce and Ljung-Box tests for autocorrelation of additional k lags with k-(p+q+1) degrees of freedom.

Forecasting

r(B)Y_t = d + q(B)e_t, t=1,2,...N

where

r(B) = 1-r₁B-r₂B²-...-r_pB^p,
q(B) = 1-q₁B-q₂B²-...-q_qB^q.

Because m = r(B)^-1d and
y(B) = r(B)^-1q(B) = 1 + y₁B + y₂B² + ...

The forecasting model can be represented as:

Y_t = m + y₁e_t-1 + y₂e_t-2 + ... + e_t

• One-Step Ahead Forecast

  Given the historical information available at the end of estimation period N,

  H_N = (Y_-p+1,...,Y_-1,Y₀,Y₁,...,Y_N; e_-q+1,...,e_-1,e₀,e₁,...,e_N)

  One-step ahead forecast is the conditional expectation of Y_N+1:

  Y_N(1) = E(Y_N+1|H_N) = m + y₁e_N + y₂e_N-1 + ...

  Compared with the observed Y_N+1 which is:

  Y_N+1 = m + e_N+1 + y₁e_N + y₂e_N-1 + ...

  One-step ahead forecast error is defined by:

  e_N(1) = Y_N+1-Y_N(1) = e_N+1

  E(e_N(1)) = 0
  s_N²(1) = Var(e_N(1)) = Var(e_N+1) = s²

• Two-Step Ahead Forecast

  Similarly, two-step ahead forecast is the conditional expectation of Y_N+2:

  Y_N(2) = E(Y_N+2|H_N) = m + y₂e_N + y₃e_N-1 + ...

  Compared with the observed Y_N+2 which is:

  Y_N+2 = m + e_N+2 + y₁e_N+1 + y₂e_N + ...

  Two-step ahead forecast error is defined by:

  e_N(2) = Y_N+2-Y_N(2) = e_N+2 + y₁e_N+1 = (1 + y₁B)e_N+2

  E(e_N(2)) = 0
  s_N²(2) = Var(e_N(2)) = Var((1+y₁B)e_N+2) = (1+y₁²)s² > s_N²(1)

• f-Step Ahead Forecast

  f-step ahead forecast is the conditional expectation of Y_N+f:

  YN(f) = E(YN+f|HN) =	m + yfeN + yf+1eN-1 + ...
  	m, as f ® ¥

  Compared with the observed Y_N+f which is:

  Y_N+f = m + e_N+f + y₁e_N+f-1 + y₂e_N+f-2 + ...

  f-step ahead forecast error is defined by:

  e_N(f) = Y_N+f-Y_N(f) = (1+y₁B+...+y_f-1B^f-1)e_N+f

  E(e_N(f)) = 0
  s_N²(f) = Var(e_N(f)) = (1+y₁²+...+y_f-1²)s² > ... > s_N²(1)

• (f+1)-Step Ahead Forecast and Forecast Revision

  In general, f+1-step ahead forecast is written as

  Y_N(f+1) = E(Y_N+f+1|H_N) = m + y_f+1e_N + y_f+2e_N-1 + ...

  Compared with the f-step ahead forecast at N+1 (with the historical information H_N+1 = (H_N,Y_N+1,e_N+1)):

  Y_N+1(f) = E(Y_N+f+1|H_N+1) = m + y_fe_N+1 + y_f+1e_N + ...

  Then the forecast revision for Y_N+f+1 is defined by

  Y_N+1(f) - Y_N(f+1) = y_fe_N+1 = y_fe_N(1)

  Therefore, with additional information available at N+1, the f-step ahead forecast is just the (f+1)-step ahead forecast made at previous date N, adjusted for one-step forecast error e_N(1) weighted by the error learning coefficient y_f as follows:

  Y_N+1(f) = Y_N(f+1) + y_fe_N(1)

Example 1

Y_t = g₀ + g₁ Y_t-1 + g₂ Y_t-2 + u_t

The alternative is to express the model (mean-deviation) residuals in terms of AR(2) structure:

Y_t = m + e_t
e_t = r₁ e_t-1 + r₂ e_t-2 + u_t

where u_t ~ nii(0,s²).

We investigate the data series Y_t for autocorrelations and partial autocorrelations. For univariate analysis, the ACF and PACF of the time series is the same as those of residuals obtained from the mean-deviation regression. Up to the maximum number of lags of ACF and PACF specified, Box-Pierece and Ljung-Box test statistics will be useful for identifying the proper order of AR(p), MA(q), or ARMA(p,q) process. In addition, Breusch-Pagan test may be used to verify the higher order of autocorrelation (Program).

Extensions

Intervention Analysis

where Z_t is a deterministic fixed variable (e.g., trend, dummy, or step variable). The interpretation of the intervention variable Z_t may be of interest.

Seasonal ARMA and Mixed Model

• Seasonal Difference

  Denote s the seasonal span (s=4 for quarterly data, s=12 for monthly data), then
  Y_t = (1-B^s)Z_t

• Seasonal ARMA

  r^s(B^s)Y_t = d + q^s(B^s)e_t

• Multiplicative Mixed Model

  r^s(B^s)r(B)Y_t = d + q^s(B^s)q(B)e_t

Transfer Function or ARMAX Model

where b(B) = b₀ + b₁B + b₂B² + ... + b_KB^K. Model analysis including model identification, estimation, and forecasting is the same as (although more complicate than) the univariate ARMA analysis. Regression parameters bs and ARMA parameters rs and qs must be simultaneously estimated through iterations of nonlinear functional (sum-of-squares or log-likelihood) optimization. For statistical reference, the degrees of freedom must be adjusted.

Multivariate Vector Autocorrelation or VAR Model