\\\\( \nonumber \newcommand{\bevisslut}{$\blacksquare$} \newenvironment{matr}[1]{\hspace{-.8mm}\begin{bmatrix}\hspace{-1mm}\begin{array}{#1}}{\end{array}\hspace{-1mm}\end{bmatrix}\hspace{-.8mm}} \newcommand{\transp}{\hspace{-.6mm}^{\top}} \newcommand{\maengde}[2]{\left\lbrace \hspace{-1mm} \begin{array}{c|c} #1 & #2 \end{array} \hspace{-1mm} \right\rbrace} \newenvironment{eqnalign}[1]{\begin{equation}\begin{array}{#1}}{\end{array}\end{equation}} \newcommand{\eqnl}{} \newcommand{\matind}[3]{{_\mathrm{#1}\mathbf{#2}_\mathrm{#3}}} \newcommand{\vekind}[2]{{_\mathrm{#1}\mathbf{#2}}} \newcommand{\jac}[2]{{\mathrm{Jacobi}_\mathbf{#1} (#2)}} \newcommand{\diver}[2]{{\mathrm{div}\mathbf{#1} (#2)}} \newcommand{\rot}[1]{{\mathbf{rot}\mathbf{(#1)}}} \newcommand{\am}{\mathrm{am}} \newcommand{\gm}{\mathrm{gm}} \newcommand{\E}{\mathrm{E}} \newcommand{\Span}{\mathrm{span}} \newcommand{\mU}{\mathbf{U}} \newcommand{\mA}{\mathbf{A}} \newcommand{\mB}{\mathbf{B}} \newcommand{\mC}{\mathbf{C}} \newcommand{\mD}{\mathbf{D}} \newcommand{\mE}{\mathbf{E}} \newcommand{\mF}{\mathbf{F}} \newcommand{\mK}{\mathbf{K}} \newcommand{\mI}{\mathbf{I}} \newcommand{\mM}{\mathbf{M}} \newcommand{\mN}{\mathbf{N}} \newcommand{\mQ}{\mathbf{Q}} \newcommand{\mT}{\mathbf{T}} \newcommand{\mV}{\mathbf{V}} \newcommand{\mW}{\mathbf{W}} \newcommand{\mX}{\mathbf{X}} \newcommand{\ma}{\mathbf{a}} \newcommand{\mb}{\mathbf{b}} \newcommand{\mc}{\mathbf{c}} \newcommand{\md}{\mathbf{d}} \newcommand{\me}{\mathbf{e}} \newcommand{\mn}{\mathbf{n}} \newcommand{\mr}{\mathbf{r}} \newcommand{\mv}{\mathbf{v}} \newcommand{\mw}{\mathbf{w}} \newcommand{\mx}{\mathbf{x}} \newcommand{\mxb}{\mathbf{x_{bet}}} \newcommand{\my}{\mathbf{y}} \newcommand{\mz}{\mathbf{z}} \newcommand{\reel}{\mathbb{R}} \newcommand{\mL}{\bm{\Lambda}} \newcommand{\mnul}{\mathbf{0}} \newcommand{\trap}[1]{\mathrm{trap}(#1)} \newcommand{\Det}{\operatorname{Det}} \newcommand{\adj}{\operatorname{adj}} \newcommand{\Ar}{\operatorname{Areal}} \newcommand{\Vol}{\operatorname{Vol}} \newcommand{\Rum}{\operatorname{Rum}} \newcommand{\diag}{\operatorname{\bf{diag}}} \newcommand{\bidiag}{\operatorname{\bf{bidiag}}} \newcommand{\spanVec}[1]{\mathrm{span}{#1}} \newcommand{\Div}{\operatorname{Div}} \newcommand{\Rot}{\operatorname{\mathbf{Rot}}} \newcommand{\Jac}{\operatorname{Jacobi}} \newcommand{\Tan}{\operatorname{Tan}} \newcommand{\Ort}{\operatorname{Ort}} \newcommand{\Flux}{\operatorname{Flux}} \newcommand{\Cmass}{\operatorname{Cm}} \newcommand{\Imom}{\operatorname{Im}} \newcommand{\Pmom}{\operatorname{Pm}} \newcommand{\IS}{\operatorname{I}} \newcommand{\IIS}{\operatorname{II}} \newcommand{\IIIS}{\operatorname{III}} \newcommand{\Le}{\operatorname{L}} \newcommand{\app}{\operatorname{app}} \newcommand{\M}{\operatorname{M}} \newcommand{\re}{\mathrm{Re}} \newcommand{\im}{\mathrm{Im}} \newcommand{\compl}{\mathbb{C}} \newcommand{\e}{\mathrm{e}} \\\\)

We continue the identification game from last time. The purpose of the exercise is to get a feeling for the iterative model building exercise that often takes place in reality when using ARIMA models. Notice that the code is changed a bit to be able to handle seasonality.

Paste the attached code into R and in groups of 2-3 do the following steps:

  1. One person describes the code line by line to the others. The rest correct the describer if he/she is wrong.
  2. Like last time, one person chooses orders for the simulated process in secret and runs the first part of the code to produce plots of the ACF and PACF of the data. Notice that there are 6 orders to choose from and a period for the seasonality. For computational reasons I suggest you to choose orders no larger than “3” are preferably smaller. Then, show the ACF and PACF to the rest of the group.
  3. The rest of the group makes a guess at an initial model structure, including all order and period of seasonality. (Notice that you might choose to only difference the data, thus choosing e.g. d=1 or D=1, while leaving the rest of the orders untouched.).
  4. The person who chose the orders for the simulation estimates a model with these orders and produdes the ACF and PACF of the residuals of this model. These are shown to the rest of the group.
  5. Iterate step 3) and 4) until a model is found for which the group is happy and evaluate whether this was close to the true model and whether you modelling approach was sensible.

You can repeat this game once or twice.

Once you are done with the game:

  1. Calculate how many different models that you would need to test, if you were to simply brute force your way through the model order identification with allowed orders up to 5.
  2. Estimate a model where p, q, P and Q are all at least 3, and keep track of roughly how much time this takes.
  3. Discuss why people often use the ACF and PACF to find model orders instead of brute force. :)
library(forecast)

p <- 1 # AR order
d <- 1 # Integration order
q <- 1 # MA order

# Seasonal part
P <- 1 # AR order
D <- 1 # Integration order
Q <- 1 # MA order
Period <- 6 # Seasonality

phi <- runif(p,-1/(p+P),1/(p+P))
theta <- runif(q,-1/(q+Q),1/(q+Q))
Phi <- runif(P,-1/(p+P),1/(p+P))
Theta <- runif(Q,-1/(q+Q),1/(q+Q))


n <- 20000
model <- Arima(ts(rnorm(n),freq=Period), order=c(p,d,q), seasonal=c(P,D,Q),
               fixed=c(phi=phi, theta=theta,
                       Phi=Phi, Theta=Theta))
Y <- simulate(model, nsim=n)

par(mfrow=c(1,2))

acf(Y,main="ACF")
pacf(Y,main="PACF")

# Choose order of estimated model

fit1 <- arima(Y,order=c(0,0,0),seasonal = list(order=c(0,0,0),period=1))

acf(residuals(fit1))
pacf(residuals(fit1))

# Iterate if necessary

cat(" Real AR parameters:",round(theta,3),"\n",
    "Estimated AR     parameters:",round(fit1<span>$coef[grep("^ar",names(fit1$</span>coef))],3),"\n",
    "Real MA parameters:",round(phi,3),"\n",
    "Estimated MA     parameters:",round(fit1<span>$coef[grep("^ma",names(fit1$</span>coef))],3),"\n",
    "Real Seasonal AR parameters:",round(theta,3),"\n",
    "Estimated Seasonal AR parameters:",round(fit1<span>$coef[grep("sar",names(fit1$</span>coef))],3),"\n",
    "Real Seasonal MA parameters:",round(phi,3),"\n",
    "Estimated Seasonal MA parameters:",round(fit1<span>$coef[grep("sma",names(fit1$</span>coef))],3))