\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{Non-stationary Modelling and Simulation} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=999) options(repos=c(CRAN="https://cloud.r-project.org/")) @ \maketitle \textbf{ARIMA Models:} An $\ARIMA(p,d,q)$ model is a model for process $\{X_t\}$ specified as \[ \Phi(\B) (1-\B)^d X_t = \Theta(\B) Z_t \] where $\Phi(\cdot)$ and $\Theta(\cdot)$ are the $p$- and $q$-order $\AR$ and $\MA$ polynomials, that is, the process $\{Y_t\}$ defined by differencing $\{X_t\}$ $d$ times at lag 1, \[ Y_t = (1-\B)^d X_t \] is a (stationary, causal) $\ARMA(p,q)$ process. Analysis of data via an $\ARIMA(p,d,q)$ model proceeds by first differencing $d$ times, and then proceeding as for an $\ARMA(p,q)$ process. Such processes are easily handled using $\texttt{R}$. For example, for the $\ARIMA(2,1,1)$ model given by the polynomials \[ \Phi(z) = (1-0.2 z)(1+0.8z) = (1+0.6z-0.16z^2) \qquad \Theta(z) = 1 + 0.9 z \] that is \[ (1+0.6 \B -0.16 \B^2) (1-\B) X_t = (1+0.9\B) Z_t \] is simulated and analyzed as follows: <>= set.seed(3834) n<-1000 X<-arima.sim(n,model=list(order=c(2,1,1),ar=c(-0.6,0.16),ma=0.9)) par(mar=c(4,4,0,0));plot(X,type='l');points(X,pch=19,cex=0.5) fit1<-arima(X,order=c(2,1,1)) fit1 @ <>= fit1.resids<-residuals(fit1) fit1.fitted<-X-fit1.resids par(mar=c(4,4,1,0));plot(X,type='l');points(X,pch=19,cex=0.5) points(fit1.fitted,pch=3,cex=0.5,col='red') legend(0,60,c('Data','Fitted'),pch=c(19,3),col=c('black','red')) round(acf(fit1.resids,lag.max=20,main='ACF for residuals',plot=F)$acf,4)[,,1] @ The residuals are uncorrelated, and plausibly white noise. The data can also be analyzed after differencing: <>= fit2<-arima(diff(X),order=c(2,0,1)) fit2 @ The parameter estimates are similar, and differ only because of the loss of a single data point that happens due to the differencing. \bigskip \textbf{Seasonal ARIMA Models:} A seasonal ARIMA (or SARIMA) model is an ARIMA model with an additional stochastic seasonal structure that can allow for seasonal differencing. The $\ARIMA(p,d,q) \times (P,D,Q)_S$ model is a model for process $\{X_t\}$ specified as \[ \Phi_S(\B^S) \Phi(\B) (1-\B^S)^D (1-\B)^d X_t = \Theta_S(\B^S) \Theta(\B) Z_t \] where the seasonal polynomials \[ \Phi_S(z) = 1 - \phi_{S1} z - \cdots - \phi_{SP} z^P \qquad \Theta_S(z) = 1 - \theta_{S1} z - \cdots - \theta_{SQ} z^Q \] allow for a seasonal ARMA structure, that is \[ \Phi_S(\B^S) = 1 - \phi_{S1} \B^S - \phi_{S2} \B^{2S} - \cdots - \phi_{SP} \B^{SP} \] and so on, as well as seasonal differencing, at lag $S$. Such models can be handled using the \texttt{sarima} package for simulation, and the \texttt{arima} function for fitting. For example, suppose we wanted to simulate and analyze the $\ARIMA(2,1,1) \times (1,0,2)_{12}$ model with \[ (1-0.3 \B^{12}) (1-1.2\B+0.8 \B^2) (1-\B) X_t = (1+0.7\B^{12}-0.3 \B^{24}) (1+0.4 \B) Z_t \] <>= library(sarima) set.seed(384) p<-2;d<-1;q<-1 P<-1;D<-0;Q<-2 S<-12 n<-1000 nvals<-n X<-sim_sarima(n=nvals, model=list(ar=c(1.2, -0.8), ma=0.4, #Non-seasonal sar=0.3, sma=c(0.7,-0.3 ), #Seasonal iorder=d, siorder=D, nseasons=S)) par(mar=c(4,4,2,0));plot(X,type='l');points(X,pch=19,cex=0.5) abline(v=100,lty=2) @ The arguments to \texttt{sim\underline{ }sarima} named \texttt{iorder} and \texttt{siorder} specify the parameters $d$ and $D$ respectively. <>= fit0<-arima(X,order=c(2,1,1),seasonal=list(order=c(1,0,2),period=12), optim.control=list(maxit=1000)) fit0 @ Note that we can combine the $\Phi_S(\cdot)$ and $\Phi(\cdot)$ filters, and the $\Theta_S(\cdot)$ and $\Theta(\cdot)$ filters, that is \begin{align*} \Phi_S(z^S) \Phi(z) & = 1 - \phi_{1}^\ast z - \cdots - \phi_{p+PS}^\ast z^{p+PS} = \Phi^\ast(z) \\[6pt] \Theta_S(z^S) \Theta(z) & = 1 - \theta_{1}^\ast z - \cdots - \theta_{q+QS}^\ast z^{q+QS} = \Theta^\ast(z) \end{align*} to obtain an $\ARMA(p+PS,q+QS)$ model before seasonal and lag one differencing. That is \[ Y_t = (1-\B^S)^D (1-\B)^d X_t \] is such that $\{Y_t\} \sim \ARMA(p+PS,q+QS)$. Accordingly, we can use the \texttt{diffinv} function to generate $\{X_t\}$ from $\{Y_t\}$ (here recall $D=0$). <>= phi.vec<-c(1.2,-0.8,rep(0,9),0.3,-0.3*1.2,0.3*0.8) theta.vec<-c(0.4,rep(0,10),0.7,0.7*0.4,rep(0,10),-0.3,-0.3*0.4) set.seed(384) Y<-arima.sim(n,model=list(order=c(p+S*P,0,q+S*Q),ar=phi.vec,ma=theta.vec)) Y<-c(diff(X[1:100]),Y);X<-diffinv(Y,lag=1);X<-X[1:n] par(mar=c(4,4,1,1));plot(X,type='l');points(X,pch=19,cex=0.5);abline(v=100,lty=2) @ Notice that simulating the data in this way cannot recover the original series without the correct `initialization'. \bigskip \textbf{Unit Roots:} In the simplest case, consider an $\AR(1)$ model \[ X_t = \phi X_{t-1} + Z_t \] from which it follows that \[ X_t - X_{t-1} = \phi X_{t-1} - X_{t-1} + Z_t = (\phi-1) X_{t-1} + Z_t \] that is \[ (1-\B) X_t = \phi^\ast X_{t-1} + Z_t \] where $\phi^\ast = \phi - 1$. This suggests an intercept-free regression relationship between the differenced data $(1-\B) X_t$ and the lagged data $X_{t-1}$. We can test the hypotheses \[ \begin{array}{lcl} H_0 & : & \phi^\ast = 0 \qquad (\text{i.e.} \: \phi = 1)\\[6pt] H_1 & : & \phi^\ast < 0 \qquad (\text{i.e.} \: \phi < 1) \end{array} \] by using the test statistic \[ t_{\text{DF}} = \frac{\widehat \phi^\ast}{\text{e.s.e}(\widehat \phi^\ast)} \] which can be computed using the \texttt{lm} function in \texttt{R}. The precise form of the test statistic is \[ \widehat \phi^\ast = \frac{\sum\limits_{t=2}^n (x_t - x_{t-1}) x_{t-1}}{\sum\limits_{t=2}^n x_{t-1}^2 } \] with estimated standard error \[ \text{e.s.e}(\widehat \phi^\ast) = \sqrt{\frac{s^2}{\sum\limits_{t=2}^n x_{t-1}^2}} \] where \[ s^2 = \frac{1}{n-2} \sum_{t=2}^n (x_t - x_{t-1} - \widehat \phi^\ast x_{t-1})^2. \] This test statistic is a form of the \textbf{Dickey-Fuller} test statistic. It can be computed using standard tools. <>= set.seed(234) n<-1000 Z<-rnorm(n) #Simulate the whote noise X<-cumsum(Z)[-n] #Create the random walk dX<-Z[-1] #Create the lagged series fit0<-lm(dX~-1+X) #Fit the regression mode coef(summary(fit0)) (phi.hat<-sum(dX*X)/sum(X^2)) #Check the calculation s<-sqrt(sum((dX-phi.hat*X)^2)/(n-2)) (ese.phi.hat<-s/sqrt(sum(X^2))) (t.df<-phi.hat/ese.phi.hat) (p.val<-pnorm(t.df)) @ The test here has non-stationarity as the null hypothesis, and the $p$-value is relatively large (i.e. \Sexpr{p.val}). However, the difficulty here is that \textit{\textbf{the test statistic is not asymptotically distributed as a standard Normal random variable}} under the null hypothesis; in fact, it is not centered at zero. This is demonstrated in the following simulation study: we simulate $N=5000$ times from a random walk and compute the test statistic each time. <>= set.seed(234) n<-1000;N<-5000 ur.test0<-function(nv){ #Function to compute the D-F test statistic. Z<-rnorm(nv);dX<-Z[-1];X<-cumsum(Z)[-n] fit0<-lm(dX~-1+X) tst<-coef(summary(fit0))[1,3] return(tst) } t.stat<-replicate(N, ur.test0(n)) #Use the replicate function for speed of computing. par(mar=c(4,4,2,0)) hist(t.stat,main='Dickey-Fuller test statistic null distribution', col='gray',xlab='t',br=seq(-5,5,by=0.25));box() #Null distribution computed via Monte Carlo abline(v=t.df,col='red',lty=2) xv<-seq(-4,4,length=1000);yv<-dnorm(xv);lines(xv,yv*0.25*N,col='red') legend(1,500,c('N(0,1)','Test statistic value'),col='red',lty=c(1,2)) @ Using this (Monte Carlo) simulation study, we estimate the $p$-value in this one-sided test by the proportion of test statistic values generated under the null hypothesis that are more negative than the observed version. <>= (p.mc<-sum(t.df > t.stat)/N) @ Thus the estimated $p$-value is \Sexpr{p.mc}. This is accurate to approximately $\sqrt{\widehat p(1-\widehat p)/N} = \Sexpr{sqrt(p.mc*(1-p.mc)/N)}$ (this is the Monte Carlo standard error), so we can conclude result is much less significant than that given initially. The empirical critical values at $0.01,0.05$ and $0.10$ are <>= quantile(t.stat,probs=c(0.01,0.05,0.1)) @ \medskip We can also use the \texttt{ur.df} function in the \texttt{urca} library: <>= library(urca) summary(ur.df(cumsum(Z),type='none',lags=0)) @ Here the critical values are computed by more extensive computation and analytical approximations. \medskip For illustration, we now consider the result for a simulated stationary process, an $\AR(1)$ model with $\phi=0.95$. <>= n<-1000 X<-arima.sim(n,model=list(ar=0.95)) dX<-diff(X) #Create the lagged series X<-X[-n] #Create the random walk fit0<-lm(dX~-1+X) #Fit the regression mode coef(summary(fit0)) (phi.hat<-sum(dX*X)/sum(X^2)) #Check the calculation s<-sqrt(sum((dX-phi.hat*X)^2)/(n-2)) (ese.phi.hat<-s/sqrt(sum(X^2))) (t.df<-phi.hat/ese.phi.hat) (p.val<-pnorm(t.df)) (p.mc<-sum(t.df > t.stat)/N) @ The test statistic is much more negative, and the Monte Carlo estimated $p$-value is now $\Sexpr{p.mc}$. Using \texttt{ur.df} we obtain <>= summary(ur.df(X,type='none',lags=0)) @ Here the test statistic is \Sexpr{t.df}, which is much more extreme than the 0.01 critical value (which is -2.58) so we reject the null hypothesis, and conclude that the process is stationary. \bigskip The Dickey-Fuller procedure can be extended to allow for a `drift' term, that is \[ X_t = \phi_0^\ast + \phi X_{t-1} + Z_t \] say, so that \[ (1-\B) X_t = \phi_0^\ast + \phi^\ast X_{t-1} + Z_t \] where we need to include an intercept in the regression model. <>= set.seed(234) library(urca) Z<-rnorm(n);X<-cumsum(Z)[-n] dX<-Z[-1] fit0<-lm(dX~X) coef(summary(fit0)) summary(ur.df(cumsum(Z),type='drift',lags=0)) @ Here we inspect the \texttt{tau2} statistic. The null distribution in this case is different again: <>= set.seed(234) n<-1000 ur.test1<-function(nv){ Z<-rnorm(nv);dX<-Z[-1];X<-cumsum(Z)[-n] fit0<-lm(dX~X) tst<-coef(summary(fit0))[2,3] return(tst) } t.stat<-replicate(N, ur.test1(n)) par(mar=c(4,4,2,2)) hist(t.stat,main='Augmented Dickey-Fuller Test Statistic', col='gray',xlab='t',br=seq(-5,5,by=0.25));box() xv<-seq(-4,4,length=1000);yv<-dnorm(xv);lines(xv,yv*0.25*N,col='red') legend(2,500,c('N(0,1)'),col='red',lty=1) @ <>= t.df<-coef(summary(fit0))[2,3] (p.mc<-sum(t.df > t.stat)/N) quantile(t.stat,probs=c(0.01,0.05,0.1)) @ \end{document}