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

\def\Ind{\mathbb{1}}



\usepackage{amsfonts,amsmath,amssymb}

\setlength{\parindent}{0pt}


\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{Spectral Estimation: Short and Long Memory Models}

\begin{document}

<<setup, cache=FALSE, echo=FALSE>>=
set.seed(3)
library(knitr)
# global chunk options
opts_chunk$set(cache=TRUE, autodep=TRUE)
options(scipen=999)
options(repos=c(CRAN="https://cloud.r-project.org/"))
myrnd <- function(x,d) {formatC(x, format="f", digits=d)}

inline_hook <- function (x) {
  if (is.numeric(x)) {
    # ifelse does a vectorized comparison
    # If integer, print without decimal; otherwise print two places
    res <- ifelse(x == round(x),
      sprintf("%6i", x),
      sprintf("%.6f", x)
    )
    paste(res, collapse = ", ")
  }
}
knit_hooks$set(inline = inline_hook)
@

\maketitle

\textbf{Model-free estimation:}  The spectral density function (SDF) can be estimated from time-domain samples.  Let $x_{1:n}$ be the time-domain data, and let
\[
\omega_k = \frac{2 \pi k}{n}
\]
for $k$ equal to
\[
-\left \lfloor \frac{n-1}{2} \right \rfloor, \ldots, \left \lfloor \frac{n}{2} \right \rfloor
\]
The set $\{\omega_k\}$ are termed the \emph{Fourier frequencies}.  For each $k$, define the $n \times 1$ vector $\bm{e}_k$ by
\[
\bm{e}_k = \frac{1}{\sqrt{n}} (e^{\imagi \omega_k}, e^{2 \imagi \omega_k},\ldots,e^{n \imagi \omega_k})^\top
\]
so that $\overline{\bm{e}}_j^\top \bm{e}_k = 0$ unless $j=k$, in which case $\overline{\bm{e}}_k^\top \bm{e}_k = 1$.  
<<se1>>=
n<-512
kvec<-seq(-floor((n-1)/2),n/2,by=1)
omega<-2*pi*kvec/n
emat<-matrix(0,nrow=n,ncol=length(omega))
for(k in 1:length(omega)){
    emat[,k]<-exp(1i*omega[k]*c(1:n))/sqrt(n)
}
Imat<-t(Conj(emat)) %*% emat
round(Imat[1:4,1:4],6)
@

$\{\bm{e}_1, \ldots, \bm{e}_n\}$ form an orthonormal basis for $\mathbb{C}^n$, that is, for any $n$-dimensional complex number $\bx$, we have
\[
\bx = \sum_{k=-\lfloor (n-1)/2 \rfloor}^{\lfloor n/2 \rfloor} a_k \bm{e}_k
\]
where
\[
a_k = \overline{\bm{e}}_k^\top \bx = \frac{1}{\sqrt{n}} \sum_{t=1}^n x_t e^{-\imagi t \omega_k}
\]
The vector $(a_1,\ldots,a_n)^\top$ is the \emph{discrete Fourier transform} of $\bx$.  Note that
\[
x_t = \sum_{k=-\lfloor (n-1)/2 \rfloor}^{\lfloor n/2 \rfloor} a_k \left[\cos (\omega_k t) + \imagi \sin(\omega_k t) \right]
\]
is the discrete version of the spectral representation theorem.

\pagebreak
\textbf{Example:} $\ARMA(1,1)$ model with $\phi=0.8$ and $\theta=0.5$.

<<se2,fig.height=5>>=
x<-arima.sim(n,model=list(ar=0.8,ma=0.5))
avec<-rep(0,n)
for(k in 1:length(omega)){
    avec[k]<-Conj(emat[,k]) %*% x
}
avec[1:5]
@

\textbf{Note:} An equivalent formulation uses the re-labelled Fourier frequencies
\[
\omega_k^\prime = \frac{2 \pi (k-1)}{n} \qquad k=1,\ldots,n
\]
and the representation using coefficients
\[
b_k = \frac{1}{\sqrt{n}} \sum_{t=1}^n x_t e^{-\imagi (t-1) \omega_k^\prime}
\]
This is the representation computed using the fast Fourier transform function \texttt{fft} in \texttt{R}, although that function does not introduce the $\sqrt{n}$ scaling.  The re-labelled Fourier frequencies are re-indexed as
\[
\omega_k^\prime = \omega_{l} \qquad l = \frac{n}{2}+k-1 \quad \mod n
\]
for $k=1,\ldots,n$.
<<se3,fig.height=5>>=
omega.new<-2*pi*c(0:(n-1))/n
bvec<-rep(0,n)
for(k in 1:length(omega)){
    bvec[k]<-sum(x*exp(-1i*c(0:(n-1))*omega.new[k]))
}
bvec[1:5]/sqrt(n)
fft(x)[1:5]/sqrt(n)
avec[(n/2)+c(0:4)]
omega[(n/2)+c(0:4)]
omega.new[1:5]
@


\pagebreak

\textbf{The Periodogram:}  The \emph{periodogram}, $\widehat I_n(\lambda)$ is defined by
\[
\widehat I_n(\lambda) = \frac{1}{n} \left| \sum_{t=1}^n x_t e^{-\imagi t \lambda } \right|^2  \qquad -\pi < \lambda \leq \pi.
\]
Note that $\widehat I_n(\omega_k) = |a_k|^2$, and
\[
\overline{\bx}^\top \bx = \sum_{k=-\lfloor (n-1)/2 \rfloor}^{\lfloor n/2 \rfloor} |a_k|^2 = \sum_{k=-\lfloor (n-1)/2 \rfloor}^{\lfloor n/2 \rfloor} \widehat I_n(\omega_k)
\]
which gives a decomposition of the total variation in $\bx$ via components at the $\omega_k$.  

<<se5>>=
sum(x^2)
sum(Mod(avec)^2)
@


Note that
\begin{eqnarray*}
\widehat I_n(\omega_k) &=& \left( \frac{1}{\sqrt{n}} \sum_{t=1}^n x_t e^{-\imagi t \omega_k} \right) \left( \frac{1}{\sqrt{n}} \sum_{s=1}^n x_s e^{\imagi s \omega_k} \right)\\[6pt]
&=& \frac{1}{n}  \sum_{t=1}^n \sum_{s=1}^n x_t x_s e^{-\imagi (t-s) \omega_k}\\[6pt]
&=& \frac{1}{n}  \sum_{t=1}^n \sum_{s=1}^n (x_t-\xmean) (x_s-\xmean) e^{-\imagi (t-s) \omega_k}\\[6pt]
&=& \sum_{h=-(n-1)}^{n-1} \widehat \gamma_n(h) e^{-\imagi h \omega_k} \qquad (h=s-t)
\end{eqnarray*}
as $\sum\limits_{t=1}^n e^{-\imagi t \omega_k} = 0$ if $\omega_k \neq 0$.  To see this last result note that
\[
\sum\limits_{t=1}^n e^{-\imagi t \omega_k}  = \frac{e^{-\imagi (n+1)\omega_k}-e^{-\imagi \omega_k} }{1-e^{-\imagi \omega_k}}
\]
by summing a finite geometric progression.  In the numerator
\begin{align*}
e^{-\imagi (n+1)\omega_k}-e^{-\imagi \omega_k} & = (\cos(2 \pi (n+1) k /n ) - \imagi \sin(2 \pi (n+1) k /n))-(\cos(2 \pi k /n) - \imagi \sin(2 \pi k /n)) \\[6pt]
& = 0.
\end{align*}
by the periodic nature $\cos$ and $\sin$.
<<se4>>=
#For example
sum(exp(-1i*c(1:n)*omega[10]))  #k=10
sum(exp(-1i*c(1:n)*omega[116])) #k=116
@


Thus $\widehat I_n(\lambda)$ is a sample-based estimate of $2 \pi f_X(\lambda)$; it is \textbf{not} consistent, however.  For $0 < \lambda_1 < \ldots < \lambda_m < \pi$,
\[
(\widehat I_n(\lambda_1),\ldots,\widehat I_n(\lambda_m))^\top \CiD (V_1,\ldots,V_n)^\top
\]
as $n \longrightarrow \infty$, where
\[
V_j \sim Exponential(1/(2 \pi f_X(\lambda_j))) .
\]
That is, as $n \longrightarrow \infty$
\[
\E[\widehat I_n(\lambda_j)] \longrightarrow 2 \pi f_X(\lambda_j)
\]
so $\widehat I_n(\lambda_j)$ is asymptotically unbiased, but $\widehat I_n(\lambda_j)$ converges in distribution, not in probability or almost surely.


\bigskip


\textbf{Example} ARMA(1,1)

\medskip

For the $\ARMA(1,1)$ model
\[
(1-0.8 B)X_t = (1+0.5B) Z_t
\]
we examine the SDF by direct computation.  We use \texttt{spectrum} to compute the periodogram: this uses the rescaled frequency scale, so we plot for ordinate $\lambda/2 \pi$.

<<sdf5, fig.height=4.5, fig.width=8>>=
library(ltsa)
acvf.vals<-tacvfARMA(phi = 0.8, theta = -0.5, maxLag = 100)
lamvec<-seq(0,0.5,by=0.001)
sdf.vec<-lamvec*0
hvec<-0:100
for(j in 1:length(lamvec)){
    sdf.vec[j]<-sum(exp(-1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[-1])+acvf.vals[1]+
                sum(exp(1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[-1])
}
par(mar=c(4,4,2,2));plot(lamvec,sdf.vec,type='l',xlab=expression(lambda/(2*pi)), 
    ylab=expression(f[X](lambda)), ylim=range(0,60),main = 'SDF')
@

It is common to plot the periodogram on the $\log_{10}$ scale, so we re-plot the SDF on this scale also.

<<sdf6, fig.height=4.25, fig.width=8>>=
par(mar=c(4,4,2,2));plot(lamvec,log10(sdf.vec),type='l',xlab=expression(lambda/(2*pi)), 
    ylab=expression(log[10](f[X](lambda))), ylim=range(-4,2),main = expression(log[10](SDF)))
@

<<sdf7, fig.height=4.25, fig.width=8>>=
set.seed(3); 
Xt<-arima.sim(n=5000,model=list(ar=0.8,ma=0.5))
par(mar=c(4,4,2,2));spectrum(Xt, main='Periodogram (log10 scale)', taper=0,detrend=FALSE)
@


<<sdf8, fig.height=4.25, fig.width=8>>=
par(mar=c(4,4,2,2));spectrum(Xt, main='Periodogram & True SDF (log10 scale)', 
    taper=0,detrend=FALSE)
lines(lamvec,sdf.vec,col='red',lwd=2)
@

For a much larger sample size, the distributional result can be more clearly seen.

<<sdf9, fig.height=4.5, fig.width=8>>=
set.seed(3); 
Xt<-arima.sim(n=50000,model=list(ar=0.8,ma=0.5))
par(mar=c(4,4,2,2));spectrum(Xt, main='Periodogram & True SDF (log10 scale, n=50000)', 
    taper=0,detrend=FALSE,xlab=expression(lambda/(2*pi)))
lines(lamvec,sdf.vec,col='red',lwd=2)
@



<<sdf10, fig.height=4.5, fig.width=8>>=
sdf.hat<-spectrum(Xt, taper=0,detrend=FALSE,plot=FALSE)$spec
n0<-length(sdf.hat)
lamvec<-c(1:n0)/(2*n0)
sdf.vec<-lamvec*0
hvec<-0:100
for(j in 1:length(lamvec)){
    sdf.vec[j]<-sum(exp(-1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[2:length(hvec)])+acvf.vals[1]+
                sum(exp(1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[2:length(hvec)])
}
par(mar=c(4,4,2,2));plot(lamvec,log10(sdf.hat),type='l',xlab=expression(lambda/(2*pi)))
lines(lamvec,log10(sdf.vec),col='red',lwd=2)
@


If in the above notation, 
\[
V_j \sim Exponential(1/(2 \pi f_X(\lambda_j))) .
\]
then
\[
\frac{V_j}{2 \pi f_X(\lambda_j)} \sim Exponential(1)
\]
which we can verify in the sample.  


<<sdf12, fig.height=4.75, fig.width=8>>=
#Compute the transformed V values
#The factor of 2 pi is already incorporated to account for the output of spectrum.
Vvec<-as.numeric(sdf.hat/sdf.vec)
mean(Vvec); var(Vvec)
par(mar=c(4,4,2,2))
hist(Vvec,nclass=100, main='Histogram of V', xlab='V', freq=FALSE)
xv<-seq(0,10,by=0.001)
lines(xv,dexp(xv),col='red',lwd=2)
legend(8,0.8,c('Density of Exp(1)'),lwd=2,col='red')
box()
@

The distributional result is maintained reasonably well for a smaller sample size also.  For $n=500$:

<<sdf13, fig.height=4.5, fig.width=8>>=
set.seed(3); 
Xt<-arima.sim(n=500,model=list(ar=0.8,ma=0.5))
sdf.hat<-spectrum(Xt, taper=0,detrend=FALSE,plot=FALSE)$spec
n0<-length(sdf.hat)
lamvec<-c(1:n0)/(2*n0)
sdf.vec<-lamvec*0
hvec<-0:100
for(j in 1:length(lamvec)){
    sdf.vec[j]<-sum(exp(-1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[2:length(hvec)])+acvf.vals[1]+
                sum(exp(1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[2:length(hvec)])
}
par(mar=c(4,4,2,2));plot(lamvec,log10(sdf.hat),type='l',
    ylab=expression(log[10](f[X](lambda))),xlab=expression(lambda/(2*pi)))
lines(lamvec,log10(sdf.vec),col='red',lwd=2)
Vvec<-as.numeric(sdf.hat/sdf.vec); c(mean(Vvec),var(Vvec))
par(mar=c(4,4,2,2)); hist(Vvec,nclass=30, main='Histogram of V', 
    xlab='V', freq=FALSE,ylim=range(0,1))
xv<-seq(0,10,by=0.001); lines(xv,dexp(xv),col='red',lwd=2)
legend(8,0.8,c('Density of Exp(1)'),lwd=2,col='red');box()
@



\textbf{Likelihood-based estimation:} For parametric models such as ARMA models with Gaussian errors, likelihood-based inference can be carried out using a large sample approximation to the time domain likelihood known as the \emph{Whittle likelihood}.  For any $\lambda$, as $n \longrightarrow \infty$
\[
\widehat I_n(\lambda) \CiD V \sim Exponential (1/(2 \pi f_X(\lambda)))
\]
If $\lambda_1,\ldots,\lambda_n$ are chosen to be the Fourier frequencies, then
\[
V_1,\ldots,V_n
\]
are approximately independent. The Whittle likelihood is then
\[
\mathcal{L}(\theta) = \prod_{j=1}^n \frac{1}{2 \pi f_X(\omega_j;\theta)} \exp\left\{ - \dfrac{\widehat I_n(\omega_j)}{2 \pi f_X(\omega_j;\theta)} \right\}.
\]
This (asymptotic) likelihood can then be maximized to estimate $\theta$.


<<sdf14, fig.height=4.5, fig.width=8>>=

log.like<-function(xp,xt){

    if(abs(xp[1])>1) return(-1e10)  #Do not allow non-stationary solution

    sdf.hat<-spectrum(xt, taper=0,detrend=FALSE,plot=FALSE)$spec
    n0<-length(sdf.hat)
    lamv<-c(1:n0)/(2*n0)
    sdf.vec<-lamv*0
    acvf.vals<-tacvfARMA(phi = xp[1], theta = -xp[2], maxLag = 100)
    hvec<-0:100
    for(j in 1:length(lamv)){
        sdf.vec[j]<-sum(exp(-1i*2*pi*lamv[j]*hvec[-1])*acvf.vals[2:length(hvec)])+
                acvf.vals[1]+
                sum(exp(1i*2*pi*lamv[j]*hvec[-1])*acvf.vals[2:length(hvec)])
    }
    Vvec<-as.numeric(sdf.hat/sdf.vec)
    return(-sum(log(sdf.vec))-sum(Vvec))
}
Xt<-arima.sim(n=5000,model=list(ar=0.8,ma=0.5))
xp0<-c(0,0)
fit.arma<-optim(xp0,fn=log.like,xt=Xt,control=list(fnscale=-1),method='BFGS')
fit.arma
@

Thus the maximizing values are $\widehat{\phi}=\Sexpr{fit.arma$par[1]}$ and $\widehat{\theta}=\Sexpr{fit.arma$par[2]}$.  

\pagebreak

For smaller sample sizes, the method is still effective: with $n=50$
<<sdf15, fig.height=4.75, fig.width=8>>=
Xt<-arima.sim(n=50,model=list(ar=0.8,ma=0.5))
xp0<-c(0,0)
fit.arma<-optim(xp0,fn=log.like,xt=Xt,control=list(fnscale=-1),method='BFGS')
fit.arma$par
fit.arma$value
acvf.vals<-tacvfARMA(phi = 0.8, theta = -0.5, maxLag = 100)
acvf.vals.fit<-tacvfARMA(phi = fit.arma$par[1], theta = -fit.arma$par[2], maxLag = 100)
lamvec<-seq(0,0.5,by=0.001)
sdf.vec<-sdf.vec.fit<-lamvec*0
hvec<-0:100
for(j in 1:length(lamvec)){
    sdf.vec[j]<-sum(exp(-1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[-1])+acvf.vals[1]+
                sum(exp(1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals[-1])
    sdf.vec.fit[j]<-sum(exp(-1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals.fit[-1])+acvf.vals.fit[1]+
                sum(exp(1i*2*pi*lamvec[j]*hvec[-1])*acvf.vals.fit[-1])
}
sdf.hat<-spectrum(Xt, taper=0,detrend=FALSE,plot=FALSE)$spec
n0<-length(sdf.hat)
lamx<-c(1:n0)/(2*n0)
par(mar=c(4,4,2,2));plot(lamx,log10(sdf.hat),type='l',ylim=range(-3,2.5),
    ylab=expression(log[10](f[X](lambda))),xlab=expression(lambda/(2*pi)))
lines(lamvec,log10(sdf.vec),col='red',lwd=2)
lines(lamvec,log10(sdf.vec.fit),col='red',lwd=2,lty=2)
legend(0.3,2.5,c('Periodogram','True','Fitted'),col=c('black','red','red'),lwd=2,lty=c(1,1,2))
@

\pagebreak

\textbf{Long memory models and Fractional Integration:} Let $\{X_t\}$ be a stationary process such that for some $\alpha \in (0,1)$ and constant $c > 0$ such that
\[
\lim_{h \longrightarrow \infty} h^\alpha \rho_X(h) =  c
\]
That is, the decay rate is \emph{polynomial} in $h$, and therefore relatively slow.  $\{X_t\}$ is a \emph{long memory} or \emph{long range dependent} process.

\bigskip

Recall the $\ARIMA(p,d,q)$ model, for $d=1,2,\ldots$,
\[
(1-B)^d \Phi(B) X_t = \Theta(B) Z_t
\]
Now suppose that $-0.5 < d < 0.5$ is a \emph{non-integer} value: note that if $d$ lies in this range, then $-d$ also lies in this range.  If $|z| < 1$, we have by the general binomial expansion, for $\alpha > -1$,

\[
(1-z)^{\alpha} = \sum_{j=0}^\infty \frac{\Gamma(j-\alpha)}{\Gamma(j+1) \Gamma(-\alpha)}  z^j
\]
where the coefficients are given by the expansion
\[
\frac{\Gamma(j-\alpha)}{\Gamma(j+1) \Gamma(-\alpha)} = \frac{(j-1-\alpha)(j-2-\alpha) \cdots (-\alpha)}{j!}  =\prod_{k=1}^j \left(\frac{k-1-\alpha}{k} \right)
\]
We apply this formula taking $\alpha=d$ for $-0.5 < d < 0.5$.

\bigskip

\textbf{Fractional White Noise:} Let $\{Z_t\} \sim \WN(0,\sigma_Z^2)$, and $W_t = (1-B)^{-d} Z_t$, so that in the original model we have
\[
\Phi(B) X_t = (1-B)^{-d} \Theta(B) Z_t = \Theta(B) (1-B)^{-d} Z_t = \Theta(B) W_t
\]
gives an $\ARMA(p,q)$ model driven by $\{W_t\}$.  $\{W_t\}$ is \emph{fractionally integrated white noise}.  In this representation, we have that for $j=1,2,\ldots$, the linear process coefficients are
\[
\psi_j = \frac{\Gamma(j+d)}{\Gamma(j+1) \Gamma(d)}
\]
$\{W_t\}$ is a stationary (linear) process.  We have that
\[
\gamma_W(h) = \sigma_Z^2 \sum_{j=0}^\infty \psi_j \psi_{j+h} = \sigma_Z^2 \sum_{j=0}^\infty \frac{\Gamma(j+d)}{\Gamma(j+1) \Gamma(d)} \frac{\Gamma(j+h+d)}{\Gamma(j+h+1) \Gamma(d)}
\] 
which can be shown (see Brockwell and Davis, \textit{Time Series: Theory and Methods}, Chapter 13) to yield
\[
\gamma_W(0) = \frac{\sigma_Z^2 \Gamma(1-2d)}{\{\Gamma(1-d)\}^2}
\]
and
\[
\rho_W(h) = \frac{\Gamma(h+d) \Gamma(1-d)}{\Gamma(h-d+1) \Gamma(d)} = \prod_{k=1}^h \left(\frac{k-1+d}{k-d}\right)
\]
We compute the infinite sum in the linear process representation for $d=0.3$ by truncating at $J=100000$ terms.
<<lm1,fig.height=4.75>>=
d<-0.3
J<-100000; jvec<-c(1:J)
psivec<-exp(lgamma(jvec+d)-lgamma(jvec+1)-lgamma(d))
psivec[1:10]; psivec[J]
H<-50; hvec<-c(1:H)
sigZsq<-2
gamma0<-sigZsq*gamma(1-2*d)/(gamma(1-d)^2)
gammah<-gamma(hvec+d)*gamma(1-d)/(gamma(hvec-d+1)*gamma(d))*gamma0
c(gamma0,gammah)


set.seed(3); n<-1000
Zt<-rnorm(n+J,0,sqrt(sigZsq))
Wt<-rep(0,n)
for(i in 1:n){
    Wt[i]<-Zt[J+i]+sum(psivec*Zt[J+i-jvec])
}
var(Wt)
par(mar=c(4,4,2,2))
plot(Wt,type='l',xlab='t', ylab=expression(W[t]))
title('Simulated realization of fractional white noise')
par(mar=c(4,4,4,2))
acf(Wt,lag.max=50,main='ACF of fractional white noise')
points(hvec,gammah/gamma0,pch=19,cex=0.7,col='red')
@

The ACVF of $\{X_t\}$ can be computed using the usual approach
\[
\gamma_X(h) = \sum_{j=0}^\infty \sum_{k=0}^\infty \psi_j \psi_k \gamma_W(h+j-k)
\]
where $\{\psi_j\}$ are the linear process coefficients corresponding to the ARMA process with polynomials $\Phi(B)$ and $\Theta(B)$.  We have by the usual result that
\[
\gamma_W(h) = \sigma_Z^2 \sum_{j=0}^\infty \psi_j \psi_{j+h}
\]

<<lm2,fig.height=4.75>>=
gammah.lp<-hvec*0
gamma0.lp<-sigZsq*(1+sum(psivec^2))
for(h in 1:H){
    ivec<-jvec+h
    ivec<-ivec[ivec <= J]
    gammah.lp[h]<-sigZsq*(psivec[h]+sum(psivec[ivec-h]*psivec[ivec]))
}
gamma0.lp
gammah.lp[1:10]
@

We can compare the computed ACVF values using the exact formula and using the truncated linear process summation.  We compare on the log scale for ease of comparison.
<<lm3,fig.height=4.75>>=
par(mar=c(4,4,2,2))
plot(c(0,hvec),log(c(gamma0,gammah)),pch=19,cex=0.5,col='red',xlab='h',ylab='log(ACVF)')
points(c(0,hvec),log(c(gamma0.lp,gammah.lp)),pch=19,cex=0.5,col='blue')
legend(30,1.0,c('Exact','Truncation'),col=c('red','blue'),pch=19)
@

These results show the limitations of studying long memory processes in the time domain: 
\begin{itemize}
  \item In the linear process representation, even out to $J = \Sexpr{J}$ terms, the coefficients are relatively large; for example, in the calculations above
  \[
  \psi_{\Sexpr{J}} = \Sexpr{psivec[J]}
  \]
  \item If we study the process using the ACVF and the sample ACVF, then we need an extremely long series to get accurate estimates.  Even with $n=\Sexpr{n}$ as in the simulation, the sample-based version is not close to the true ACVF.
       
\end{itemize}

However, in the frequency domain, analysis is a little more straightforward.  By direct calculation using the linear process representation and the general result
\[
f_X(\lambda) = \frac{\sigma_Z^2}{2 \pi} |\Psi(e^{-\imagi \lambda}) |^2
\]
we have here that
\[
f_W(\lambda) = \frac{\sigma_Z^2}{2 \pi} \frac{1}{|(1-e^{-\imagi \lambda})^{d} |^2} = \frac{\sigma_Z^2}{2 \pi} \frac{1}{|1-e^{-\imagi \lambda}|^{2d}}
\]

<<lm4,fig.height=4.75>>=
n<-500
Zt<-rnorm(n+J,0,sqrt(sigZsq))
Wt<-rep(0,n)
for(i in 1:n){
    Wt[i]<-Zt[J+i]+sum(psivec*Zt[J+i-jvec])
}
var(Wt)
sdf.hat<-spectrum(Wt, taper=0,detrend=FALSE,plot=FALSE)$spec
lamvec<-seq(0,0.5,by=0.001)
sdf.vec<-sigZsq*Mod(1-exp(-1i*2*pi*lamvec))^{-2*d}
n0<-length(sdf.hat)
lamx<-c(1:n0)/(2*n0)
par(mar=c(4,4,2,2));plot(lamx,log10(sdf.hat),type='l',ylim=range(-3,2.5),
    ylab=expression(log[10](f[X](lambda))),xlab=expression(lambda/(2*pi)))
lines(lamvec,log10(sdf.vec),col='red',lwd=2)
@



We can again estimate parameters $d$ and $\sigma_Z^2$ using the Whittle likelihood.

<<lm5, fig.height=4.5, fig.width=8>>=

log.like.fwn<-function(xp,xt){

    if(abs(xp[1])>0.5) return(-1e10)  #Do not allow non-stationary solution
    dv<-xp[1]
    if(xp[2]<=0) return(-1e10)
    sdf.hat<-spectrum(xt, taper=0,detrend=FALSE,plot=FALSE)$spec
    n0<-length(sdf.hat)
    lamv<-c(1:n0)/(2*n0)
    sdf.vec<-xp[2]*Mod(1-exp(-1i*2*pi*lamv))^{-2*dv}
    Vvec<-as.numeric(sdf.hat/sdf.vec)
    return(-sum(log(sdf.vec))-sum(Vvec))
}
xp0<-c(0.3,2)
fit.fwn<-optim(xp0,fn=log.like.fwn,xt=Wt,control=list(fnscale=-1),method='Nelder-Mead')
fit.fwn$par
fit.fwn$value
@
Hence the estimate of $d$ is \Sexpr{fit.fwn$par[1]}, and the estimate of $\sigma_Z^2$ is $\Sexpr{fit.fwn$par[2]}$.

\medskip

\textbf{Limiting behaviour:} If $\{X_t\} \sim \ARIMA(p,d,q)$ with $-1/2 < d < 1/2$, then
\[
f_X(\lambda) = \frac{\sigma_Z^2}{2 \pi} \frac{\left|\Theta(e^{-\imagi \lambda})\right|^2}{\left|\Phi(e^{-\imagi \lambda})\right|^2} \frac{1}{|1-e^{-\imagi \lambda}|^{2d}}
\]
where, if $p=q=0$,
\[
f_X(\lambda) = f_W(\lambda) = \frac{\sigma_Z^2}{2 \pi} \frac{1}{|1-e^{-\imagi \lambda}|^{2d}}
\]
In this case note that, by a Taylor-type expansion
\[
|1-e^{-\imagi \lambda}|^{2d} = \bigg| 1 - \left(1-\imagi \lambda + \frac{\lambda^2}{2} + \text{o}(\lambda^2) \right) \bigg|^{2d} = |\lambda+\text{o}(\lambda)|^{2d}
\]
as $\lambda \longrightarrow 0$.  Thus for $\lambda$ near zero
\[
f_X(\lambda) \approx \frac{\sigma_Z^2}{2 \pi |\lambda|^{2d}}
\]
and the SDF is unbounded at zero.  If $\{W_t\} \sim \ARIMA(0,d,0)$, that is,
\[
W_t = (1-B)^{-d} Z_t
\]
then using the binomial expansion method
\begin{itemize}

\item $\MA(\infty)$ representation
\[
W_t = \sum_{j=0}^\infty \psi_j Z_{t-j} \qquad \psi_j = \frac{\Gamma(j+d)}{\Gamma(j+1) \Gamma(d)}
\]

\item $\AR(\infty)$ representation
\[
Z_t = \sum_{j=0}^\infty \pi_j W_{t-j} \qquad \pi_j = \frac{\Gamma(j-d)}{\Gamma(j+1) \Gamma(-d)}
\]

\end{itemize}
For large $j$
\begin{itemize}

\item $\MA(\infty)$ representation
\[
\psi_j \approx \frac{j^{d-1}}{\Gamma(d)}
\]

\item $\AR(\infty)$ representation
\[
\pi_j \approx \frac{j^{-d-1}}{\Gamma(-d)}
\]

\end{itemize}



If $\{X_t\} \sim \ARIMA(p,d,q)$, write
\[
(1-B)^{d} \Phi(B) X_t = \Theta(B) Z_t
\]
and suppose $0 < d < 1/2$. Then for large $h$
\[
\gamma_X(h) \approx C |h|^{2d - 1} \propto \frac{1}{|h|^{1-2d}}
\]
The implication is that for $h>0$ large enough,
\[
\log \gamma_X(h) \approx \beta_0 - (1-2 d ) \log h
\]
so $\log \gamma_X(h)$ is approximately linear in $\log h$, with slope $2d-1$.  Similarly, for the $p=q=0$ case
\[
\frac{\widehat{I}_j(\omega_j)}{2 \pi f_X(\lambda)} = \frac{\widehat{I}_j(\omega_j) |1-e^{-\imagi \lambda}|^{2d}}{\sigma_Z^2} = U_j \sim Exponential(1)
\]
say. Therefore
\[
\log \widehat{I}_j(\omega_j) + 2d \log |1-e^{-\imagi \lambda}| - \log \sigma_Z^2 = \log U_j
\]
and hence
\[
\log \widehat{I}_j(\omega_j) = \log \sigma_Z^2 - 2d \log |1-e^{-\imagi \lambda}| + \log U_j
\]
Making the approximation $|1-e^{-\imagi \lambda}| \approx |\lambda|$, we have that for $\lambda$ near zero
\[
\log \widehat{I}_j(\omega_j) = \log \sigma_Z^2 - 2d \log |\lambda| + \log U_j
\]
so plotting the log periodogram against log $\lambda$ should yield a straight line with slope $-2d$.

<<lm6, fig.height=4.5, fig.width=8>>=
n0<-length(sdf.hat)
lamx<-c(1:n0)/(2*n0)
par(mar=c(4,4,2,2));plot(log(lamx[1:50]),log(sdf.hat[1:50]),ylim=range(-5,5),
     pch=19,cex=0.6,ylab=expression(log(f[X](lambda))),xlab=expression(log(lambda/(2*pi))))
x<-log(lamx[1:50])
y<-log(sdf.hat[1:50])
abline(lm(y~x)$coef,col='red',lty=2)
summary(lm(y~x))$coef
@

For a larger sample size
<<lm7, fig.height=4.5, fig.width=8>>=
n<-1000
Zt<-rnorm(n+J,0,sqrt(sigZsq))
Wt<-rep(0,n)
for(i in 1:n){
    Wt[i]<-Zt[J+i]+sum(psivec*Zt[J+i-jvec])
}
sdf.hat<-spectrum(Wt, taper=0,detrend=FALSE,plot=FALSE)$spec
n0<-length(sdf.hat)
lamx<-c(1:n0)/(2*n0)

par(mar=c(4,4,2,2));plot(log(lamx[1:100]),log(sdf.hat[1:100]),ylim=range(-5,5),
    pch=19,cex=0.6,ylab=expression(log(f[X](lambda))),xlab=expression(log(lambda/(2*pi))))
x<-log(lamx[1:100])
y<-log(sdf.hat[1:100])
abline(lm(y~x)$coef,col='red',lty=2)
summary(lm(y~x))$coef
@
\end{document}