\documentclass[notitlepage]{article}
\usepackage{../../Math545}
\usepackage{listings}
\usepackage{numprint}
\usepackage{enumitem}
\usepackage{bbm}

\usepackage{amsfonts,amsmath}

%\def\E{\Expect}

\setlength\parindent{0in}

\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{Decomposition}

\begin{document}

\maketitle

We seek to decompose an observed series $\{Y_t\}$ as
\[
Y_t = m_t + s_t + X_t
\]
where
\begin{itemize}

\item $m_t$ is a deterministic trend

\item $s_t$ is a seasonal component with period $d$

\item $X_t$ is a zero mean stationary process.
\end{itemize}

\bigskip

\textbf{Smoothing:}  For a non-seasonal data series $\{x_t\}$, we may proceed by smoothing to remove non-linear trends:

\begin{itemize}

\item \textbf{Smoothing by averaging:}  Suppose
\[
\widehat m_t = \frac{1}{2q+1} \sum_{i=-q}^q x_{t+i}
\]
with $x_s = 1$ for $s < 1$ and $s >T$.  More generally,
\[
\widehat m_t = \sum_{i=-q}^q a_i x_{t+i}
\]
for constants $a_{-q},\ldots,a_q$.



\item \textbf{Exponential Smoothing:}  Suppose $\widehat m_1 = x_1$, and for $t \geq
2$
\[
\widehat m_t = \alpha x_t + (1-\alpha) \widehat m_{t-1}
\]
for $0 < \alpha < 1$. It follows that
\[
\widehat m_t = (1-\alpha)^{t-1} x_1 + \sum_{j=1}^{t-2} \alpha (1-\alpha)^j x_{t-j} + \alpha x_t
\]
If $\alpha$ is near one, there is strong dependence on recent $x$ values.

\end{itemize}

\bigskip

\textbf{EXAMPLE:} Suppose we have the process
\[
X_t = m_t + E_t
\]
for
\[
m_t = 10 + 0.01 t + 5 \exp\{-0.001(t-100)^2\}
\]
and $E_t \sim \AR(1)$ with parameter $\phi_1 = 0.6$ based on a Gaussian IID process with variance $\sigma_Z^2 = 2^2$

<<d1,fig.height=5,fig.width=8>>=
set.seed(43)
n<-500
tv<-seq(1,n,by=1)
true.m<-10+0.01*tv+5*exp(-0.0001*(tv-100)^2)
phi<-0.6
sigZ<-2
et<-xt<-rep(0,n)
et[1]<-rnorm(1)*sigZ/sqrt(1-phi^2)
xt[1]<-true.m[1]+et[1]
for(i in 2:n){
 et[i]<-phi*et[i-1]+rnorm(1)*sigZ
 xt[i]<-true.m[i]+et[i]   
}
plot(tv,xt,type='l',xlab='t',ylab=expression(x[t]),ylim=range(5,25))
points(tv,xt,pch=19,cex=0.5)
lines(tv,true.m,col='red',lwd=2)
legend(350,25,c('True mean'),lty=1,lwd=2,col='red')
@

We attempt smoothing by averaging with $q=5$, and exponential smoothing with $\alpha=0.1$.

<<d2,fig.height=4.75,fig.width=8>>=
m1<-m2<-rep(0,n)

#Averaging
q<-5
for(i in (q+1):(n-q)){
    window<-(i-q):(i+q)
    m1[i]<-mean(xt[window])
}

#Exponential smoothing
alpha<-0.1
m2[1]<-xt[1]
for(i in 2:n){
    m2[i]<-alpha*xt[i]+(1-alpha)*m2[i-1]
}
par(mar=c(4,4,2,2))
plot(tv,xt,type='l',xlab='t',ylab=expression(x[t]),ylim=range(5,25))
points(tv,xt,pch=19,cex=0.5)
lines(tv[(q+1):(n-q)],m1[(q+1):(n-q)],col='blue',lwd=2)
legend(350,25,c('Averaging'),lty=1,lwd=2,col='blue')

plot(tv,xt,type='l',xlab='t',ylab=expression(x[t]),ylim=range(5,25))
points(tv,xt,pch=19,cex=0.5)
lines(tv,m2,col='green',lwd=2)
legend(300,25,c('Exponential smoothing'),lty=1,lwd=2,col='green')
@

\textbf{General Strategy:} For series $\{Y_t\}$
\[
Y_t = m_t + s_t + X_t
\]
where $s_t$ has period $d$.

\begin{enumerate}[label=\arabic*.]

\item If $d$ is even, set $q=d/2$, else set $q = (d-1)/2$.

\item Use moving average for detrending respecting the seasonality
\[
\widehat m_t = \left\{
\begin{array}{ll}
\dfrac{y_{t-q}+y_{t+q}}{2d} + \dfrac{1}{d} \sum\limits_{j=-q+1}^{q-1} y_{t+j} & d\text{ even}\\
\\
\dfrac{1}{d} \sum\limits_{j=-q}^{q} y_{t+j} & d\text{ odd}\\
\end{array}
\right.
\]



\item Remove seasonality: for $k=1,\ldots,d$, compute
\[
w_k = \frac{1}{n_d} \sum_{j} (y_{k+jd}-\widehat m_{k+jd})
\]
where $n_d$ is the number of terms in the sum (essentially the number of
cycles in the series).  For $k=1,\ldots,d$, set
\[
\widehat s_k = w_k - \frac{1}{d}\sum_{j=1}^d w_j
\]
and set $\widehat s_k = \widehat s_{k-d}$ for $k=d+1,d+2,\ldots,T$.  The deseasonalized data are $y_t^\star = y_t - \widehat s_t$.



\item Recompute the deseasonalized trend $\widehat m_t^\star$ from the deasonalized data $y_t^\star$ using the procedure in Step 2.

\item Compute the residual series
\[
\widehat x_t = y_t - \widehat m_t^\star - \widehat s_t.
\]
\end{enumerate}

The resulting process $\widehat x_t$ should be a zero mean, but potentially autocorrelated process.

\bigskip

\textbf{EXAMPLE:} Suppose we have the process
\[
Y_t = m_t + s_t + X_t
\]
for
\[
m_t = 10 + 0.01 t + 5 \exp\{-0.001(t-100)^2\}
\]
and 
\[
s_t = 2 \cos(2 \pi \lambda t)
\]
so that $d=1/\lambda$.  In the simulation below, we take $\lambda = 1/52$, so $d=52$.  We again set $X_t \sim \AR(1)$ with parameter $\phi_1 = 0.6$ based on a Gaussian IID process with variance $\sigma_Z^2 = 2^2$

<<d3,fig.height=5,fig.width=8>>=
set.seed(43)
n<-500
tv<-seq(1,n,by=1)
true.m<-10+0.01*tv+5*exp(-0.0001*(tv-100)^2)
lambda<-1/52
d<-1/lambda
true.s<-2*cos(2*pi*lambda*tv)
phi<-0.6
sigZ<-2
yt<-xt<-rep(0,n)
xt[1]<-rnorm(1)*sigZ/sqrt(1-phi^2)
yt[1]<-true.m[1]+true.s[1]+xt[1]
for(i in 2:n){
 xt[i]<-phi*xt[i-1]+rnorm(1)*sigZ
 yt[i]<-true.m[i]+true.s[i]+xt[i]   
}
plot(tv,yt,type='l',xlab='t',ylab=expression(x[t]),ylim=range(5,25))
points(tv,yt,pch=19,cex=0.5)
lines(tv,true.m+true.s,col='red',lwd=2)
legend(350,25,c('True mean'),lty=1,lwd=2,col='red')
@

We implement the strategy as follows:

\begin{enumerate}[label=\arabic*.]

\item Here $d=52$, so we set $q=d/2 = 26$.

\item For each $t$, we use the averaging formula
\[
\widehat m_t = \dfrac{y_{t-q}+y_{t+q}}{2d} + \dfrac{1}{d} \sum\limits_{j=-q+1}^{q-1} y_{t+j}
\]


<<d4,fig.height=5,fig.width=8>>=
q<-d/2
mt<-rep(0,n)
for(i in 1:q){
    window<-(i-q+1):(i+q-1)
    window<-window[window>0]
    mt[i]<-mean(yt[window])
}
for(i in (q+1):(n-q)){
    mt[i]<-sum(yt[(i-q+1):(i+q-1)])/d+(yt[i-q]+yt[i+q])/(2*d)
}
for(i in (n-q+1):n){
    window<-(i-q+1):(i+q-1)
    window<-window[window<=n]
    mt[i]<-mean(yt[window])
}
plot(tv,mt,type='l',xlab='t',lwd=2,ylab=expression(m[t]),ylim=range(5,25),col='red')
points(tv,yt,pch=19,cex=0.5)
legend(300,25,c('Initial mean estimate'),lty=1,lwd=2,col='red')
@


\item For $k=1,\ldots,d$, compute
\[
w_k = \frac{1}{n_d} \sum_{j} (y_{k+jd}-\widehat m_{k+jd})
\]
where $n_d$ is the number of terms in the sum (essentially the number of cycles in the series).  For $k=1,\ldots,d$, set
\[
\widehat s_k = w_k - \frac{1}{d}\sum_{j=1}^d w_j
\]
and set $\widehat s_k = \widehat s_{k-d}$ for $k=d+1,d+2,\ldots,T$.  The deseasonalized data are $y_t^\star = y_t - \widehat s_t$.

<<d5,fig.height=5,fig.width=8>>=
wk<-rep(0,d)
for(k in 1:d){
    indices<-k+c(1:11)*d
    indices<-indices[indices<=n]
    wk[k]<-mean(yt[indices]-mt[indices])
}
sk<-wk-mean(wk)
par(mar=c(4,5,2,2))
plot(1:d,sk,xlab='k',pch=19,cex=0.5,ylab=expression(s[k]),ylim=range(-5,5),col='blue')
legend(25,5,c('Estimated seasonal terms'),pch=19,col='blue')
st<-rep(sk,11)
st<-st[1:n]
plot(tv,st,xlab='t',pch=19,cex=0.5,type='l',ylab=expression(s[t]),ylim=range(-5,5),col='blue')
legend(300,5,c('Estimated seasonal component'),pch=19,col='blue')
ystar<-yt-st
plot(tv,ystar,type='l',xlab='t',lwd=2,ylab=expression(paste(y[t]^"*")),ylim=range(5,25))
points(tv,ystar,pch=19,cex=0.5)
legend(300,25,c('Deseasonalized data'),lty=1,lwd=2,pch=19)
@


\item Recompute the deseasonalized trend $\widehat m_t^\star$ from the deasonalized data $y_t^\star$ using the procedure in Step 2.

<<d6,fig.height=5,fig.width=8>>=
mstar<-rep(0,n)
for(i in 1:q){
    window<-(i-q+1):(i+q-1)
    window<-window[window>0]
    mstar[i]<-mean(ystar[window])
}
for(i in (q+1):(n-q)){
    mstar[i]<-sum(ystar[(i-q+1):(i+q-1)])/d+(ystar[i-q]+ystar[i+q])/(2*d)
}
for(i in (n-q+1):n){
    window<-(i-q+1):(i+q-1)
    window<-window[window<=n]
    mstar[i]<-mean(yt[window])
}
plot(tv,mstar,type='l',xlab='t',lwd=2,ylab=expression(paste(m[t]^"*")),
    ylim=range(5,25),col='red')
points(tv,ystar,pch=19,cex=0.5)
legend(300,25,c('Final mean estimate'),lty=1,lwd=2,col='red')
@

\item Compute the residual series
\[
\widehat x_t = y_t - \widehat m_t^\star - \widehat s_t.
\]

We can examine the properties of this residual series
<<d7,fig.height=4.5,fig.width=8>>=
xf<-yt-mstar-st
par(mar=c(4,4,1,2))
plot(tv,xf,type='l',xlab='t',ylab=expression(x[t]),ylim=range(-10,10))
points(tv,xf,pch=19,cex=0.5); 
legend(200,10,c('Final de-trended/de-seasonalized data'),lty=1)
@

The residual series seems to be zero mean and constant variance, and the estimated ACF sequence resembles that of a $\AR(1)$ process, although there is still some apparent excess correlation around lag $d=52$.
<<d8,fig.height=5,fig.width=8>>=
par(mar=c(4,4,4,2))
acf(xf,lag.max=60,main='ACF of residual series')
@


\end{enumerate}

\bigskip
\textbf{Real data example:}  Air Passenger data

These data represent the counts of the number (in thousands) of air passengers in the US between January 1949 and December 1960. 

<<d9,fig.height=5.0,fig.width=8>>=
Air<-read.csv('https://www.math.mcgill.ca/dstephens/545/Data/airline-passengers.csv')
head(Air,n=4L)
AirPassengers<-Air[,2]
y<-ts(AirPassengers,start=c(1949,1),frequency=12)
n<-length(y)
par(mar=c(4,4,2,2))
plot(y,type="l",ylim=range(0,700),ylab="Number of Passengers (thousands)",xlab="Months")
lines(y)
points(y,pch=19,cex=0.5)
legend(0,600,c("Data"),lty=c(1,1),col=c("black"),pch=c(19))
@

The data are recorded monthly, and there is an apparent annual cycle.  We follow the procedure above, with $d=12$ (and hence $\lambda = 1/12$.

<<d10,fig.height=5,fig.width=8>>=
yt<-as.numeric(Air[,2])
tv<-1:n
d<-12
q<-d/2
mt<-rep(0,n)
for(i in 1:q){
    window<-(i-q+1):(i+q-1)
    window<-window[window>0]
    mt[i]<-mean(yt[window])
}
for(i in (q+1):(n-q)){
    mt[i]<-sum(yt[(i-q+1):(i+q-1)])/d+(yt[i-q]+yt[i+q])/(2*d)
}
for(i in (n-q+1):n){
    window<-(i-q+1):(i+q-1)
    window<-window[window<=n]
    mt[i]<-mean(yt[window])
}
par(mar=c(4,4,2,2))
plot(tv,mt,type='l',xlab='t',lwd=2,ylab=expression(m[t]),ylim=range(0,700),col='red')
points(tv,yt,pch=19,cex=0.5)
legend(0,700,c('Initial mean estimate'),lty=1,lwd=2,col='red')

wk<-rep(0,d)
for(k in 1:d){
    indices<-k+c(1:21)*d
    indices<-indices[indices<=n]
    wk[k]<-mean(yt[indices]-mt[indices])
}
sk<-wk-mean(wk)
sk #Estimated seasonal components

st<-rep(sk,21)
st<-st[1:n]; ystar<-yt-st
plot(tv,st,type='l',xlab='t',pch=19,cex=0.5,ylab=expression(s[t]),ylim=range(-100,100),col='blue')
points(tv,st,pch=19,cex=0.5,col='blue')
legend(0,100,c('Estimated seasonal component'),pch=19,col='blue')
plot(tv,ystar,type='l',xlab='t',pch=19,cex=0.5,ylab=expression(paste(y[t]^"*")),ylim=range(0,700))
points(tv,ystar,pch=19,cex=0.5)
legend(0,700,c('Deseasonalized data'),pch=19)
mstar<-rep(0,n)
for(i in 1:q){
    window<-(i-q+1):(i+q-1);  window<-window[window>0]
    mstar[i]<-mean(ystar[window])
}
for(i in (q+1):(n-q)){
    mstar[i]<-sum(ystar[(i-q+1):(i+q-1)])/d+(ystar[i-q]+ystar[i+q])/(2*d)
}
for(i in (n-q+1):n){
    window<-(i-q+1):(i+q-1);     window<-window[window<=n]
    mstar[i]<-mean(yt[window])
}
plot(tv,mstar,type='l',xlab='t',lwd=2,ylab=expression(paste(m[t]^"*")),
        ylim=range(0,700),col='red')
points(tv,ystar,pch=19,cex=0.5)
legend(0,700,c('Final mean estimate'),lty=1,lwd=2,col='red')
xf<-yt-mstar-st
par(mar=c(2,4,3,2))
plot(tv,xf,type='l',xlab='t',ylab=expression(x[t]),ylim=range(-100,100))
points(tv,xf,pch=19,cex=0.5)
legend(0,100,c('Final de-trended/de-seasonalized data'),lty=1)
acf(xf,lag.max=60,main='ACF of residual series')
@

\pagebreak

We can see here that the seasonality has not been removed successfully.  This is because the amplitude of the seasonal variation is not constant.  This can be partly addressed by taking a log transform of the data.

<<d11,fig.height=5,fig.width=8>>=
y<-ts(log(AirPassengers),start=c(1949,1),frequency=12)
n<-length(y)
par(mar=c(4,4,2,2))
plot(y,type="l",ylim=range(4,7),ylab="Number of Passengers (thousands) - log scale",xlab="Months")
lines(y)
points(y,pch=19,cex=0.5)
legend(0,7,c("Data"),lty=c(1,1),col=c("black"),pch=c(19))
@

<<d12,fig.height=5,fig.width=8>>=
yt<-as.numeric(log(Air[,2]))
tv<-1:n
d<-12
q<-d/2
mt<-rep(0,n)
for(i in 1:q){
    window<-(i-q+1):(i+q-1)
    window<-window[window>0]
    mt[i]<-mean(yt[window])
}
for(i in (q+1):(n-q)){
    mt[i]<-sum(yt[(i-q+1):(i+q-1)])/d+(yt[i-q]+yt[i+q])/(2*d)
}
for(i in (n-q+1):n){
    window<-(i-q+1):(i+q-1)
    window<-window[window<=n]
    mt[i]<-mean(yt[window])
}
par(mar=c(4,4,2,2))
plot(tv,mt,type='l',xlab='t',lwd=2,ylab=expression(m[t]),ylim=range(4,7),col='red')
points(tv,yt,pch=19,cex=0.5)
legend(0,7,c('Initial mean estimate'),lty=1,lwd=2,col='red')

wk<-rep(0,d)
for(k in 1:d){
    indices<-k+c(1:21)*d
    indices<-indices[indices<=n]
    wk[k]<-mean(yt[indices]-mt[indices])
}
sk<-wk-mean(wk)
sk #Estimated seasonal components

st<-rep(sk,21)
st<-st[1:n]
plot(tv,st,type='l',xlab='t',pch=19,cex=0.5,ylab=expression(s[t]),ylim=range(-1,1),col='blue')
points(tv,st,pch=19,cex=0.5,col='blue')
legend(0,1,c('Estimated seasonal component'),pch=19,col='blue')
ystar<-yt-st
plot(tv,ystar,type='l',xlab='t',pch=19,cex=0.5,ylab=expression(paste(y[t]^"*")),ylim=range(4,7))
points(tv,ystar,pch=19,cex=0.5); legend(0,7,c('Deseasonalized data'),pch=19)
mstar<-rep(0,n)
for(i in 1:q){
    window<-(i-q+1):(i+q-1);  window<-window[window>0]
    mstar[i]<-mean(ystar[window])
}
for(i in (q+1):(n-q)){
    mstar[i]<-sum(ystar[(i-q+1):(i+q-1)])/d+(ystar[i-q]+ystar[i+q])/(2*d)
}
for(i in (n-q+1):n){
    window<-(i-q+1):(i+q-1);     window<-window[window<=n]
    mstar[i]<-mean(yt[window])
}
plot(tv,mstar,type='l',xlab='t',lwd=2,ylab=expression(paste(m[t]^"*")),
        ylim=range(4,7),col='red')
points(tv,ystar,pch=19,cex=0.5)
legend(0,7,c('Final mean estimate'),lty=1,lwd=2,col='red')
xf<-yt-mstar-st
par(mar=c(2,4,3,2))
plot(tv,xf,type='l',xlab='t',ylab=expression(x[t]),ylim=range(-0.15,0.15))
points(tv,xf,pch=19,cex=0.5)
legend(0,0.15,c('Final de-trended/de-seasonalized data'),lty=1)
acf(xf,lag.max=60,main='ACF of residual series')
@

We observe more plausibly stationary residual series, although there is some indication of non-constant variance.
\end{document}