Book Notes: Intro to Time Series and Forecasting -- Ch6 ARIMA Models

For the pdf slides, click here

When data is not stationary

  • Implication of not stationary: sample ACF or sample PACF do not rapidly decrease to zero as lag increases

  • What shall we do?
    • Differencing, then fit an ARMA →→ ARIMA
    • Transformation, then fit an ARMA
    • Seasonal model →→ SARIMA

A non-stationary exmaple: Dow Jones utilities index data

library(itsmr); ## Load the ITSM-R package
par(mfrow = c(1, 3));
plot.ts(dowj, main = 'Raw data'); 
acf(dowj); pacf(dowj);

After differencing

par(mfrow = c(1, 3));
dowj_diff = dowj[-length(dowj)] - dowj[-1];
plot.ts(dowj_diff, main = 'Data after differencing'); 
acf(dowj_diff); pacf(dowj_diff);

ARIMA Models

ARIMA model: definition

  • Autoregressive integrated moving-average models (ARIMA): Let d∈Nd∈N, then series {Xt}{Xt} is an ARIMA(p,d,q)(p,d,q) process if Yt=(1−B)dXtYt=(1−B)dXt is a causal ARMA(p,q)(p,q) process.

  • Difference equation (DE) for an ARIMA(p,d,q)(p,d,q) process ϕ∗(B)Xt=ϕ(B)(1−B)dXt=θ(B)Zt,{Zt}∼WN(0,σ2)ϕ∗(B)Xt=ϕ(B)(1−B)dXt=θ(B)Zt,{Zt}∼WN(0,σ2)
    • ϕ(z)ϕ(z): polynomial of degree pp, and ϕ(z)≠0ϕ(z)≠0 for |z|≤1|z|≤1
    • θ(z)θ(z): polynomial of degree qq
    • ϕ∗(z)=ϕ(z)(1−z)dϕ∗(z)=ϕ(z)(1−z)d: has a zero of order dd at z=1z=1

  • An ARIMA process with d>0d>0 is NOT stationary!

ARIMA mean is not dertermined by the DE

  • {Xt}{Xt} is an ARIMA(p,d,q)(p,d,q) process. We can add an arbitrary polynomial trend of degree d−1d−1 Wt=Xt+A0+A1t+⋯+Ad−1td−1Wt=Xt+A0+A1t+⋯+Ad−1td−1 with A0,…,Ad−1A0,…,Ad−1 being any random variables, and {Wt}{Wt} still satisfies the same ARIMA(p,d,q)(p,d,q) difference equation

  • In other words, the ARIMA DE determines the second-order properties of {(1−B)dXt}{(1−B)dXt} but not those of {Xt}{Xt}

    • For parameter estimation: ϕ, θ, and σ2 are estimated based on {(1−B)dXt} rather than {Xt}
    • For forecast, we need additional assumptions

Fit data using ARIMA processes

  • Whether to fit a finite time series using
    • non-stationary models (such as ARIMA), or
    • directly using stationary models (such as ARMA)?
  • If the fitted stationary ARMA model’s ϕ(⋅) have zeros very close to unit circles, then fitting an ARIMA model is better
    • Parameter estimation is stable
    • The differenced series may only need a low-order ARMA
  • Limitation of ARIMA: only permits data to be nonstationary in a very special way
    • Non-stationary: can have zeros anywhere on the unit circle |z|=1
    • ARIMA model: only has a zero of multiplicity d at the point z=1

Transformation and Identification Techniques

Natural log transformation

  • When data variance increases with mean, it’s common to apply log transformation before fitting the data using ARIMA or ARMA.

  • When does log transfomation work well? Suppose that E(Xt)=μt,Var(Xt)=σ2μ2t Then by first-order Taylor expansion of log(Xt) at μt: log(Xt)≈log(μt)+Xt−μtμt ⟹ Var[log(Xt)]≈Var(Xt)μ2t=σ2 The data after log transformation log(Xt) has a constant variance

  • Note: log transformation can only be applied to positive data

  • Note: If Yt=log(Xt), then because expectation and logarithm are not interchangeable, ˆXt≠exp(ˆYt)

Generalize the log transformation: Box-Cox transformation

  • Box-Cox transformation fλ(x)={xλ−1λ,x≥0,λ>0log(x),x>0,λ=0

    • Usual range: 0≤λ≤1.5
    • Common values: λ=0,0.5

  • Note: limλ→0fλ(x)=log(x)

  • Box-Cox transformation can only be applied to non-negative data

Unit Root Test

Unit root test for AR(1) process

  • {Xt} is an AR(1) process: Xt−μ=ϕ1(Xt−1−μ)+Zt

  • Equivalent DE: ∇Xt=Xt−Xt−1=ϕ∗0+ϕ∗1Xt−1+Zt
    • where ϕ∗0=μ(1−ϕ1) and ϕ∗1=ϕ1−1
    • Regressing ∇Xt onto 1 and Xt−1, we get the OLS estimator ˆϕ∗1 and its standard error SE(ˆϕ∗1)
  • Augmented Dickey-Fuller test for AR(1)
    • Hypotheses: H0:ϕ1=1 ⟷ H1:ϕ1<1
    • Equivalent hypotheses: H0:ϕ∗1=0 ⟷ H1:ϕ∗1<0
    • Test statistic: limit distribution under H0 is not normal or t ˆτ=ˆϕ∗1SE(ˆϕ∗1)
    • Rejection region: reject H0 if {ˆτ<−2.57(90%)ˆτ<−2.86(95%)ˆτ<−3.43(99%)

Unit root test for AR(p) process

  • AR(p) process: Xt−μ=ϕ1(Xt−1−μ)+⋯+ϕp(Xt−p−μ)+Zt

  • Equivalent DE: ∇Xt=ϕ∗0+ϕ∗1Xt−1+ϕ∗2∇Xt−1+⋯+ϕ∗p∇Xt−p+1+Zt

    • where ϕ∗0=μ(1−∑pi=1ϕi), ϕ∗1=∑pi=1ϕi−1, and ϕ∗j=−∑pi=jϕi for j≥2
    • Regressing ∇Xt onto 1,Xt−1,∇Xt−1,⋯,∇Xt−p+1, we get the OLS estimator ˆϕ∗1 and its standard error SE(ˆϕ∗1)
  • Augmented Dickey-Fuller test for AR(p)
    • Hypotheses: H0:ϕ∗1=0 ⟷ H1:ϕ∗1<0
    • Test statistic: ˆτ=ˆϕ∗1SE(ˆϕ∗1)
    • Rejection region: same as augmented Dickey-Fuller test for AR(1)

Implement augmented Dickey-Fuller test in R

library(tseries);
## Note: the lag k here is the AR order p
adf.test(dowj, k = 2);
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dowj
## Dickey-Fuller = -1.3788, Lag order = 2, p-value = 0.8295
## alternative hypothesis: stationary

Forecast ARIMA models

Forecast an ARIMA(p,1,q) process

  • {Xt} is ARIMA(p,1,q), and {Yt=∇Xt} is a causal ARMA(p,q) Xt=X0+t∑j=1Yj,t=1,2,…

  • Best linear predictor of Xn+1 PnXn+1=Pn(X0+Y1+⋯+Yn+1)=Pn(Xn+Yn+1)=Xn+Pn(Yn+1),
    • Pn means based on {1,X0,X1,…,Xn}, or equivalently, {1,X0,Y1,…,Yn}
    • To find Pn(Yn+1), we need to know E(X20) and E(X0Yj), for j=1,…,n+1.

  • A sufficient assumption for Pn(Yn+1) to be the best linear predictor in term of {Y1,…,Yn}: X0 is uncorrelated with Y1,Y2,…

Forecast an ARIMA(p,d,q) process

  • The observed ARIMA(p,d,q) process {Xt} satisfies Yt=(1−B)dXt,t=1,2,…,{Yt}∼causal ARMA(p,q)

  • Assumption: the random vector (X1−d,…,X0) is uncorrelated with Yt for all t>0

  • One-step predictors ˆYn+1=PnYn+1 and ˆXn+1=PnXn+1: Xn+1−ˆXn+1=Yn+1−ˆYn+1

  • Recall: the h-step predictor of ARMA(p,q) for n>max(p,q): PnYn+h=p∑i=1ϕiPnYn+h−i+q∑j=hθn+h−1,j(Yn+h−j−ˆYn+h−j)

  • h-step predictor of ARIMA(p,d,q) for n>max(p,q): PnXn+h=p+d∑i=1ϕ∗iPnXn+h−i+q∑j=hθn+h−1,j(Xn+h−j−ˆXn+h−j) where ϕ∗(z)=(1−z)dϕ(z)=1−ϕ∗1z−⋯−ϕ∗p+dzp+d

Seasonal ARIMA Models

Seasonal ARIMA (SARIMA) Model: definition

  • Suppose d,D are non-negative integers. {Xt} is a seasonal ARIMA(p,d,q) × (P,D,Q)s process with period s if the differenced series Yt=(1−B)d(1−Bs)DXt is a causal ARMA process defined by ϕ(B)Φ(Bs)Yt=θ(B)Θ(Bs)Zt,{Zt}∼WN(0,σ2) where ϕ(z)=1−ϕ1z−⋯−ϕpzp,Φ(z)=1−Φ1z−⋯−ΦPzPθ(z)=1+θz+⋯+θqzq,Θ(z)=1+Θz+⋯+ΘQzQ

  • {Yt} is causal if and only if neither ϕ(z) or Φ(z) has zeros inside the unit circle

  • Usually, s=12 for monthly data

Special case: seasonal ARMA (SARMA)

  • Between-year model: for monthly data, each one of the 12 time series is generated by the same ARMA(P,Q) model Φ(B12)Yt=Θ(B12)Ut,{Uj+12t,t∈Z}∼WN(0,σ2U)

  • SARMA(P,Q) with period s: in the above between-year model, the period 12 can be changed to any positive integer s
    • If {Ut,t∈Z}∼WN(0,σ2U), then the ACVF γ(h)=0 unless h divides s evenly. But this may not be ideal for real life application! E.g., this Feb is correlated with last Feb, but not this Jan.
  • General SARMA(p,q) × (P,Q) with period s: incorporate dependency between the 12 series by letting {Ut} be ARMA: Φ(Bs)Yt=Θ(Bs)Ut,ϕ(B)Ut=θ(B)Zt,{Zt}∼WN(0,σ2)
    • Equivalent DE for the general SARMA: ϕ(B)Φ(Bs)Yt=θ(B)Θ(Bs)Zt

Fit a SARIMA Model

  • Period s is known

  1. Find d and D to make the difference series {Yt} to look stationary

  2. Examine the sample ACF and sample PACF of {Yt} at lags being multiples of s, to find orders P,Q

  3. Use ˆρ(1),…,ˆρ(s−1) to find orders p,q

  4. Use AICC to decide among competing order choices

  5. Given orders (p,d,q,P,D,Q), estimate MLE of parameters (ϕ,θ,Φ,Θ,σ2)

Regression with ARMA Errors

Regression with ARMA errors: OLS estimation

  • Linear model with ARMA errors W=(W1,…,Wn)′: Yt=x′tβ+Wt,t=1,…,n,{Wt}∼causal ARMA(p,q)
    • Note: each row is indexed by a different time t!
    • Error covariance Γn=E(WW′)
  • Ordinary least squares (OLS) estimator ˆβOLS=(X′X)−1X′Y
    • Estimated by minimizing (Y−Xβ)′(Y−Xβ)
    • OLS is unbiased, even when errors are dependent!

Regression with ARMA errors: GLS estimation

  • Generalized least squares (GLS) estimator ˆβGLS=(X′Γ−1nX)−1X′Γ−1nY

    • Estimated by minimizing the weighted sum of squares (Y−Xβ)′Γ−1n(Y−Xβ)

    • Covariance Cov(ˆβGLS)=(X′Γ−1nX)−1

    • GLS is the best linear unbiased estimator, i.e., for any vector c and any unbiased estimator ˆβ, we have Var(c′ˆβGLS)≤Var(c′ˆβ)

When {Wt} is an AR(p) process

  • We can apply ϕ(B) to each side of the regression equation and get uncorrelated, zero-mean, constant-variance errors ϕ(B)Y=ϕ(B)Xβ+ϕ(B)W=ϕ(B)Xβ+Z

  • Under the transformed target variable Y∗t=ϕ(B)Yt,t=p+1,…,n and the transformed design matrix X∗t,j=ϕ(B)Xt,j,t=p+1,…,n,j=1,…,k the OLS estimator is the best linear unbiased estimator

  • Note: after the transformation, the regression sample size reduces to n−p

Regression with ARMA errors: MLE

  • MLE of β,ϕ,θ,σ2 can be estimated by maximizing the Gaussian likelihood with error covariance Γn

  • An iterative scheme

    1. Compute ˆβOLS and regression residuals Y−XˆβOLS

    2. Based on the estimated residuals, compute MLE of the ARMA(p,q) parameters

    3. Based on the fitted ARMA model, compute ˆβGLS

    4. Compute regression residuals Y−XˆβGLS, and return to Step 2 until estimators stabilize

  • Asymptotic properties of MLE: If {Wt} is a causal and invertible ARMA, then
    • MLEs are asymptotically normal
    • Estimated regression coefficients are asymptotically independent of estimated ARMA parameters

References