\documentclass[notitlepage]{article} \usepackage{../../Math545} \usepackage{listings} \usepackage{numprint} \usepackage{enumerate} \usepackage{bbm} \usepackage{bm} \usepackage{amsfonts,amsmath,amssymb} \newtheorem{alg}{Algorithm} \newtheorem{ex}{Example} \newtheorem{note}{Note} \newenvironment{proof}{\par\noindent{\normalfont{\bf{Proof }}}} {\hfill \small$\blacksquare$\\[2mm]} \setlength{\parindent}{0in} \def\E{\mathbb{E}} \def\Var{\text{Var}} \def\Cov{\text{Cov}} \def\Corr{\text{Corr}} \def\bW{\mathbf{W}} \newcommand{\bs}[1]{\bm{#1}} \newcommand{\Rho}{\mathrm{P}} \newcommand{\cif}{\text{ if }} \newcommand{\Ra}{\Longrightarrow} \newcommand{\ra}{\longrightarrow} \def\bphi{\boldsymbol \phi} \lstloadlanguages{R} \definecolor{keywordcolor}{rgb}{0,0.6,0.6} \definecolor{delimcolor}{rgb}{0.461,0.039,0.102} \definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432} \lstdefinestyle{Rsettings}{ basicstyle=\ttfamily, breaklines=true, showstringspaces=false, keywords={if, else, function, theFunction, tmp}, % Write as many keywords otherkeywords={}, commentstyle=\itshape\color{Rcommentcolor}, keywordstyle=\color{keywordcolor}, moredelim=[s][\color{delimcolor}]{"}{"}, } \lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}} \def\ctitle{MATH 545: Introduction to Time Series Analysis} \def\ttitle{Winter 2025} \def\stitle{Selecting a Time Series Model} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) @ \maketitle <>= set.seed(2384) n<-30*12 Y<-as.numeric(arima.sim(n,model=list(ar=c(0.5,-0.56)),rand.gen = function(n, ...) sqrt(0.5) * rnorm(n))) Y<-cumsum(Y) @ We consider a systematic approach to analyzing time series data using the following example. Here is a basic summary of the data. <>= length(Y) Y[1:20] par(mar=c(4,4,1,0)) plot(Y,pch=19,cex=0.5);lines(Y) abline(h=mean(Y),lty=2) @ There is no obvious systematic trend or periodicity, as indicated by the deviations about the plotted mean. <>= mean(Y) @ However, there does seem to be strong positive dependence in the series. This can be assessed by inspecting the sample acf: <>= par(mar=c(4,4,1,0)) acf(Y,lag.max=50) @ This plot (with a slow decay of the acf) indicates that a form of non-stationarity may be present. Recall that there are different ways that non-stationarity can be present: \begin{itemize} \item deterministic variation in mean or covariance; \item deterministic periodic or seasonal component; \item stochastic non-stationarity due to a \textit{unit root} in an $\AR$ component (i.e. $|\alpha_j| = 1$), leading to random walk-like behaviour. \end{itemize} In this case, the last of these seems most likely, so we attempt to remove it using first differencing. <>= X<-diff(Y) par(mar=c(4,4,1,0)) plot(X,pch=19,cex=0.5);lines(X);abline(h=0,lty=2) acf(X,ylim=range(-1,1)) #ACF pacf(X,ylim=range(-1,1)) #PACF @ These plots indicate that after differencing, the process seems stationary; the acf plot indicates that the acf diminishes to zero rapidly, and the pacf plot indicates that the pacf is only non-zero for a small number of lags. Therefore an $\AR(p)$ model (possibly an $\ARMA(p,q)$ model) with $p$ (and $q$) small should be adequate. \medskip We can then use model selection using AIC or BIC to assess which exact model might be best. Recall that the AIC and BIC calculations essentially assume that the data are Gaussian, so that the log-likelihood for data $\bx_{1:n} = (x_1,\ldots,x_n)^\top$ is \[ \ell_n(\theta,\phi,\sigma_Z) = \text{const.} - \frac{1}{2} \log |\Gammat_X(n;\theta,\phi,\sigma_Z) | - \frac{1}{2} \bx_{1:n}^\top \{\Gammat_X(n;\theta,\phi,\sigma_Z) \}^{-1} \bx_{1:n}. \] We then perform estimation by maximizing this function with respect to the parameters to obtain estimates $(\widehat \theta,\widehat \phi, \widehat \sigma_Z)$. The AIC and BIC are then defined by \begin{itemize} \item AIC: \[ \text{AIC}_{n}(p,q) = - 2 \ell_n(\widehat \theta,\widehat \phi, \widehat \sigma_Z) + 2 (p+q+1) \] \item BIC: \[ \text{BIC}_{n}(p,q) = - 2 \ell_n(\widehat \theta,\widehat \phi, \widehat \sigma_Z) + (\log n) (p+q+1) \]\end{itemize} We select the model for which the criterion is \textbf{minimized}. First we compare $\AR(p)$ models for $p=0,1,2,\ldots,10$: <>= fit.0<-ar(X,order.max=10,demean=T,method='mle') fit.0 fit.0$aic @ The AIC values reported here are the values above the minimizing value. These results indicate that an $\AR(5)$ model is selected. We can attempt the same analysis using the \texttt{arima} function <>= AIC.vec<-BIC.vec<-rep(0,11) for(p in 0:10){ fit.p<-arima(X,order=c(p,0,0),include.mean=T,method='ML') AIC.vec[1+p]<-fit.p$aic BIC.vec[1+p]<-BIC(fit.p) } cbind(0:p,AIC.vec,BIC.vec) par(mar=c(4,4,1,0)) plot(0:p,AIC.vec,pch=19,col='blue',ylim=range(c(AIC.vec,BIC.vec)),xlab='p',ylab='AIC or BIC') lines(0:p,AIC.vec,col='blue') points(0:p,BIC.vec,pch=19,col='red');lines(0:p,BIC.vec,col='red') legend(8,900,c("AIC","BIC"),pch=19,lty=1,col=c('blue','red')) abline(v=which.min(AIC.vec)-1,col='blue',lty=2) abline(v=which.min(BIC.vec)-1,col='red',lty=2) @ Here the AIC and BIC values are on an absolute scale (that is, not relative to the minimizing value). The AIC is minimized at $p=5$ as before, but the BIC is minimized at $p=2$. Overall, there seems to be little to choose between the $\AR(p)$ models for $p=2,3,4,5$. \medskip We can also perform the calculation for $\ARMA(p,q)$ models. <>= AIC.mat<-BIC.mat<-matrix(0,7,7) for(p in 0:6){ for(q in 0:6){ fit.pq<-arima(X,order=c(p,0,q),include.mean=T,method='ML',optim.control = list(maxit=1000)) AIC.mat[1+p,1+q]<-fit.pq$aic BIC.mat[1+p,1+q]<-BIC(fit.pq) } } rownames(AIC.mat)<-paste("p=",c(0:6),sep="");colnames(AIC.mat)<-paste("q=",c(0:6),sep="") rownames(BIC.mat)<-paste("p=",c(0:6),sep="");colnames(BIC.mat)<-paste("q=",c(0:6),sep="") @ Again, we can inspect the AIC and BIC values relative to the minimum value: <>= round(AIC.mat-min(AIC.mat),4) @ For the AIC there are several competing models: the minimum is obtained for the $\ARMA(2,5)$ model, however, there are seven other models within 1.00 of this minimizing value. We can compare the fitted acfs of the following models <>= fit.25<-arima(X,order=c(2,0,5),include.mean=T,method='ML',optim.control = list(maxit=1000)) fit.53<-arima(X,order=c(5,0,3),include.mean=T,method='ML',optim.control = list(maxit=1000)) fit.56<-arima(X,order=c(5,0,6),include.mean=T,method='ML',optim.control = list(maxit=1000)) AIC(fit.25,fit.53,fit.56) acf.25<-ARMAacf(ar=coef(fit.25)[1:2],ma=coef(fit.25)[3:7],lag.max=30) acf.53<-ARMAacf(ar=coef(fit.53)[1:5],ma=coef(fit.53)[6:8],lag.max=30) acf.56<-ARMAacf(ar=coef(fit.56)[1:5],ma=coef(fit.56)[6:11],lag.max=30) par(mar=c(4,4,3,0)) acf(X,lag.max=30,ylim=range(-1,1),main='Sample and fitted ACFs') points(0:30,acf.25,pch=19,cex=0.6,col='blue') points(0:30,acf.53,pch=19,cex=0.6,col='red') points(0:30,acf.56,pch=19,cex=0.6,col='green') legend(20,1,c('ARMA(2,5)','ARMA(5,3)','ARMA(5,6)'),col=c('blue','red','green'),pch=19) @ The fitted ACF values are very similar, and match the sample acf values reasonably well. \medskip For the BIC, the selected model is the $\MA(3)$, and only the $\ARMA(2,1)$ model is close in BIC terms. <>= round(BIC.mat-min(BIC.mat),4) @ We may compute the residual series from the fits of these models and compare them <>= fit.03<-arima(X,order=c(0,0,3),include.mean=T,method='ML',optim.control = list(maxit=1000)) fit.21<-arima(X,order=c(2,0,1),include.mean=T,method='ML',optim.control = list(maxit=1000)) res.03<-residuals(fit.03);res.21<-residuals(fit.21) par(mar=c(4,4,3,0)) acf(res.03,lag.max=30,ylim=range(-1,1),main='ACF of MA(3) residuals') acf(res.21,lag.max=30,ylim=range(-1,1),main='ACF of ARMA(2,1) residuals') @ Both of the residual series seem to be plausibly white noise processes. \medskip Finally, we might compare forecasts from several of these models: for this, we can use the \texttt{predict} function: <>= pred.25<-predict(fit.25,n.ahead=10)$pred pred.03<-predict(fit.03,n.ahead=10)$pred pred.21<-predict(fit.21,n.ahead=10)$pred par(mar=c(6,4,1,0)) plot(340:359,X[340:359],xlim=range(340,370),pch=19,cex=0.6,xlab='Time',ylab='X') lines(340:359,X[340:359]) points(360:369,pred.25,col='blue',pch=19,cex=0.6);lines(359:369,c(X[359],pred.25),col='blue') points(360:369,pred.03,col='red',pch=19,cex=0.6);lines(359:369,c(X[359],pred.03),col='red') points(360:369,pred.21,col='green',pch=19,cex=0.6);lines(359:369,c(X[359],pred.21),col='green') abline(h=0,lty=2) legend(360,1.5,c('ARMA(2,5)','MA(3)','ARMA(2,1)'),col=c('blue','red','green'),pch=19,lty=1) @ The forecasts are very similar overall due to the fact that the autocovariances for the different models are similar. As usual, the forecasts return to the process mean of zero as $h$ increases. \end{document}