\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 and frequency domain concepts} \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/")) 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 The ACVF of a (zero-mean) stationary process $\{X_t\}$ \[ \gamma_X(h) = \E[X_t X_{t+h}] \] is a symmetric, non-negative definite function of the discrete argument $h = 0, \pm 1, \pm 2, \ldots$. We exploit the fact that $\gamma_X(h)$ has a unique representation in terms of \emph{frequency contributions}. Let $\{X_t\}$ be stationary and zero-mean, with ACVF $\gamma_X(.)$. Then the \emph{spectral density function} (SDF) of $\{X_t\}$ is the function $f_X(.)$ defined by \[ f_X(\lambda) = \frac{1}{2 \pi} \sum_{h=-\infty}^\infty e^{-\imagi \lambda h} \gamma_X(h) \qquad \lambda \in \R \] provided this sum is convergent; if \[ \sum_{h=-\infty}^\infty |\gamma_X(h)| < \infty. \] then $f_X(\lambda)$ is well-defined, as $|e^{-\imagi \lambda h}| = 1$, so \[ \sum_{h=-\infty}^\infty |e^{-\imagi \lambda h} \gamma_X(h)| \leq \sum_{h=-\infty}^\infty |e^{-\imagi \lambda h}|| \gamma_X(h)| = \sum_{h=-\infty}^\infty |\gamma_X(h)| < \infty. \] although the absolute summability is not a necessary condition. Also, $f_X(\lambda)$ is \emph{periodic} as \[ e^{-\imagi \lambda h} = \cos(\lambda h) - \imagi \sin (\lambda h) \] so we can restrict attention to $\lambda \in (-\pi,\pi]$. We have that \begin{itemize} \item[(i)] $f_X(\lambda) = f_X(-\lambda)$ and {\small \begin{align*} f_X(\lambda) = \frac{1}{2 \pi} \sum_{h=-\infty}^\infty e^{-\imagi \lambda h} \gamma_X(h) & = \frac{1}{2 \pi} \left[ \gamma_X(0) + \sum_{h=1}^\infty (e^{-\imagi \lambda h} + e^{\imagi \lambda h}) \gamma_X(h) \right] & \text{as } \gamma(-h) = \gamma(h)\\[6pt] & = \frac{1}{2 \pi} \left[ \gamma_X(0) + \sum_{h=1}^\infty 2 \cos(\lambda h) \gamma_X(h) \right] \\[6pt] & = \frac{1}{2 \pi} \sum_{h=-\infty}^\infty \cos(\lambda h) \gamma_X(h) & \text{as } \cos(-\lambda h) = \cos(\lambda h) \end{align*}} \item[(ii)] $f_X(\lambda) \geq 0$, $-\pi < \lambda \leq \pi$. \item[(iii)] Inversion formula: \begin{align*} \int_{-\pi}^\pi e^{\imagi h \lambda} f_X(\lambda) \: d \lambda & = \frac{1}{2 \pi} \int_{-\pi}^\pi \left\{ \sum_{k=-\infty}^\infty e^{ \imagi \lambda (h-k)} \gamma_X(k) \: \right\} d \lambda= \frac{1}{2 \pi} \sum_{k=-\infty}^\infty \gamma_X(k) \left\{\int_{-\pi}^\pi e^{ \imagi \lambda (h-k)} \: d \lambda \right\} = \gamma_X(h). \end{align*} \end{itemize} The last result follows as if $k \neq h$, setting $j = h-k$ \begin{align*} \int_{-\pi}^\pi e^{ \imagi \lambda (h-k)} \: d \lambda & = \int_{-\pi}^\pi [\cos(\lambda j) + \imagi \sin (\lambda j)] \: d \lambda = \int_{-\pi}^\pi \cos(\lambda j) \: d \lambda + \imagi \int_{-\pi}^\pi \sin (\lambda j) \: d \lambda = 0. \end{align*} Note also that \[ \gamma_X(h) = \int_{-\pi}^\pi e^{\imagi h \lambda} f_X(\lambda) \: d \lambda = \int_{-\pi}^\pi \left[\cos(h \lambda) + \imagi \sin (h \lambda) \right] f_X(\lambda) \: d \lambda = \int_{-\pi}^\pi \cos(h \lambda) f_X(\lambda) \: d \lambda \] as $\sin(x)$ is an odd function. \pagebreak There is a 1-1 correspondence between ACVFs and spectral density functions; if \begin{itemize} \item[(i)] $f_X(\lambda) \geq 0$, for all $\lambda \in (0,\pi]$. \item[(ii)] For $h=0,\pm 1,\pm 2,\ldots$, if \[ \gamma_X(h) = \int_{-\pi}^\pi e^{\imagi h \lambda} f_X(\lambda) \: d \lambda \] \end{itemize} then $f_X$ is the SDF of $\{X_t\}$ with ACVF $\gamma_X(h)$. For a generalization, a function $\gamma(.)$, $\gamma : \mathbb{Z} \longrightarrow \R$, is an autocovariance function if and only if there exists a right continuous, non-decreasing bounded function $F$ on $[-\pi,\pi]$ with $F(-\pi) = 0$, and for all $h \in \mathbb{Z}$ \[ \gamma (h) = \int_{(-\pi,\pi]} e^{\imagi h \lambda} \: dF(\lambda). \] As the corresponding series $\{\gamma(h)\}$ is real-valued, $F$ satisfies \[ \int_{(a,b]} dF(x) = \int_{(-b,-a]} dF(x) \] \begin{itemize} \item The notation $dF(\lambda)$ represents the infinitesimal change in $F$ at $\lambda$. \item $F$ is a \emph{generalized distribution function} as $F(\lambda)/F(\pi)$ is a distribution function on $[-\pi,\pi]$. \item $F$ is termed the \emph{spectral distribution function}; it can be purely discrete, absolutely continuous, or mixed. If $F$ is absolutely continuous, then the \emph{spectral density function} $f$ is defined by \[ F(\lambda) = \int_{-\pi}^\lambda f(x) \: dx \qquad \lambda \in (-\pi,\pi] \] \end{itemize} \bigskip \textbf{Harmonic Process:} Suppose \[ X_t = U \cos (\omega t) + V \sin (\omega t) \] where $U$ and $V$ are uncorrelated, with \[ \E[U] = \E[V] = 0 \qquad \qquad \Var[U] = \Var[V] = 1 \] Then by a previous result \[ \gamma_X(h) = \cos (\omega h) \qquad h \in \mathbb{Z} \] so that \[ \sum_{h=-\infty}^\infty |\gamma_X(h)| = \sum_{h=-\infty}^\infty |\cos (\omega h)| = \infty \] and $|\gamma_X(h)|$ does not converge to any limit. However, let \[ F_X(\lambda) = \left\{ \begin{array}{lc} 0 & \lambda < -\omega\\[6pt] \dfrac{1}{2} & -\omega \leq \lambda < \omega\\[9pt] 1 & \omega \leq \lambda \end{array} \right. \] and define \[ f_X(\lambda) = \left\{ \begin{array}{lc} \dfrac{1}{2} & \lambda = -\omega\\[9pt] \dfrac{1}{2} & \lambda = \omega\\[9pt] 0 & \text{otherwise} \end{array} \right. \] with the interpretation \[ \int_{(-\pi,\pi)} e^{\imagi h \lambda} \: dF_X(\lambda) = \frac{1}{2} e^{-\imagi h \lambda} + \frac{1}{2} e^{\imagi h \lambda} = \cos(h \omega) \] Thus $F_X(\lambda)$ is a right continuous step function with steps at $ \pm \omega$, and $f_X(\lambda)$ takes the value of $1/2$ at $ \pm \omega$. The variation in $X$ has a characteristic frequency $\omega$ that indicates its periodicity. <>= lam<-seq(-pi,pi,by=0.01) omega<-pi/2 Flam<-(lam >= -omega)*0.5+(lam >= omega)*0.5 Flam[which.min((lam+omega)^2)]<-NA Flam[which.min((lam-omega)^2)]<-NA par(mar=c(4,4,2,2)) plot(lam, Flam,type='l',xlab=expression(lambda),ylab=expression(F[X](lambda))) points(lam[which.min((lam+omega)^2)+1],Flam[which.min((lam+omega)^2)+1],pch=19) points(lam[which.min((lam-omega)^2)-1],Flam[which.min((lam-omega)^2)+1],pch=19) plot(lam,lam*0,type='n',xlab=expression(lambda),ylab=expression(f[X](lambda)),ylim=range(0,1)) lines(rep(lam[which.min((lam+omega)^2)],2),c(0,1/2)); points(lam[which.min((lam+omega)^2)],1/2,pch=19) lines(rep(lam[which.min((lam-omega)^2)],2),c(0,1/2)); points(lam[which.min((lam-omega)^2)],1/2,pch=19) @ \textbf{Superposition of harmonic processes:} Suppose \[ X_t = \sum_{j=1}^k \left[ U_j \cos (\omega_j t) + V_j \sin (\omega_j t) \right] \] where $U_1,\ldots,U_k$ and $V_1,\ldots, V_k$ are uncorrelated, with \[ \E[U_j] = \E[V_j] = 0 \qquad \qquad \Var[U_j] = \Var[V_j] = \sigma_j^2 \] and where $0 < \omega_1 < \cdots < \omega_k < \pi$. Then \[ \gamma_X(h) = \sum_{j=1}^k \sigma_j^2 \cos(\omega_j h) \qquad h \in \mathbb{Z}. \] Then \[ F_X(\lambda) = \sum_{j=1}^k \sigma_j^2 F_j(\lambda) \] where \[ F_j(\lambda) = \left\{ \begin{array}{lc} 0 & \lambda < -\omega_j\\[6pt] \dfrac{1}{2} & -\omega_j \leq \lambda < \omega_j\\[12pt] 1 & \omega_j \leq \lambda \end{array} \right. \] Note that \[ \Var[X_t] = \gamma_X(0) = \sum_{j=1}^k \sigma_j^2 \] Suppose $k=3$, with \begin{itemize} \item $j=1$: $\omega_1 = \pi/4$, $\sigma_1^2 = 0.4$ \item $j=2$: $\omega_2 = \pi/3$, $\sigma_2^2 = 1.2$ \item $j=3$: $\omega_3 = 2\pi/3$, $\sigma_3^2 = 0.7$ \end{itemize} <>= lam<-seq(-pi,pi,by=0.01) omega<-c(pi/4,pi/3, 2*pi/3) sigsq<-c(0.4,1.2,0.7) hvec<-c(0:25) gammah<-hvec*0 for(j in 1:3){ gammah<-gammah+sigsq[j]*cos(omega[j]*hvec) } acvf<-cbind(hvec,gammah) rownames(acvf)<-paste('h=',hvec,sep='') colnames(acvf)<-c('h','ACVF') acvf[1:10,2] @ <>= Flam<-lam*0 par(mar=c(4,4,2,2)) cvec<-c('red','green','blue','black') plot(lam, Flam,type='n',xlab=expression(lambda),ylab=expression(F[X](lambda)),ylim=range(0,sum(sigsq))) for(j in 1:3){ Flamj<-sigsq[j]*((lam >= -omega[j])*0.5+(lam >= omega[j])*0.5) Flam<-Flam+Flamj Flamj[which.min((lam+omega[j])^2)]<-NA Flamj[which.min((lam-omega[j])^2)]<-NA Flam[which.min((lam+omega[j])^2)]<-NA Flam[which.min((lam-omega[j])^2)]<-NA lines(lam, Flamj,type='l',col=cvec[j]) points(lam[which.min((lam+omega[j])^2)+1],Flamj[which.min((lam+omega[j])^2)+1],pch=19,col=cvec[j]) points(lam[which.min((lam-omega[j])^2)-1],Flamj[which.min((lam-omega[j])^2)+1],pch=19,col=cvec[j]) } legend(-3,2.0,c('j=1','j=2','j=3','Sum'),lty=1,col=cvec) lines(lam,Flam) for(j in 1:3){ points(lam[which.min((lam+omega[j])^2)+1],Flam[which.min((lam+omega[j])^2)+1],pch=19) points(lam[which.min((lam-omega[j])^2)-1],Flam[which.min((lam-omega[j])^2)+1],pch=19) } plot(lam,lam*0,type='n',xlab=expression(lambda),ylab=expression(f[X](lambda)),ylim=range(0,sum(sigsq))) for(j in 1:3){ lines(rep(lam[which.min((lam+omega[j])^2)],2),c(0,sigsq[j])); points(lam[which.min((lam+omega[j])^2)],sigsq[j],pch=19) lines(rep(lam[which.min((lam-omega[j])^2)],2),c(0,sigsq[j])); points(lam[which.min((lam-omega[j])^2)],sigsq[j],pch=19) } par(mar=c(4,4,2,2)) plot(hvec,gammah,type='p',pch=19,xlab='h',ylab='ACVF',ylim=range(-3,3)) for(j in 1:length(hvec)){lines(c(hvec[j],hvec[j]),c(0,gammah[j]))} k<-3 Uvec<-rnorm(k)*sqrt(sigsq) Vvec<-rnorm(k)*sqrt(sigsq) tvec<-seq(0,500,by=1) Xvec<-tvec*0 for(j in 1:k){ Xvec<-Xvec+Uvec[j]*cos(omega[j]*tvec)+Vvec[j]*sin(omega[j]*tvec) } par(mar=c(4,4,2,2)) plot(tvec,Xvec,type='l',xlab='t',ylab=expression(X[t])) @ \pagebreak Now suppose $k=50$, with \begin{itemize} \item $\omega_j = j \pi/k$, \item $\sigma_j^2 = j/k$ \end{itemize} for $j=1,\ldots,k$. <>= k<-50 lam<-seq(-pi,pi,by=0.01) omega<-c(1:k)*pi/k sigsq<-c(1:k)/k hvec<-c(0:25) gammah<-hvec*0 for(j in 1:k){ gammah<-gammah+sigsq[j]*cos(omega[j]*hvec) } acvf<-cbind(hvec,gammah) rownames(acvf)<-paste('h=',hvec,sep='') colnames(acvf)<-c('h','ACVF') acvf[1:10,2] @ <>= Flam<-lam*0 par(mar=c(4,4,2,2)) cvec<-c('red','green','blue','black') plot(lam, Flam,type='n',xlab=expression(lambda),ylab=expression(F[X](lambda)),ylim=range(0,sum(sigsq))) for(j in 1:k){ Flamj<-sigsq[j]*((lam >= -omega[j])*0.5+(lam >= omega[j])*0.5) Flam<-Flam+Flamj Flamj[which.min((lam+omega[j])^2)]<-NA Flamj[which.min((lam-omega[j])^2)]<-NA Flam[which.min((lam+omega[j])^2)]<-NA Flam[which.min((lam-omega[j])^2)]<-NA } lines(lam,Flam) for(j in 1:k){ points(lam[which.min((lam+omega[j])^2)+1],Flam[which.min((lam+omega[j])^2)+1],pch=19,cex=0.5) points(lam[which.min((lam-omega[j])^2)-1],Flam[which.min((lam-omega[j])^2)+1],pch=19,cex=0.5) } plot(lam,lam*0,type='n',xlab=expression(lambda),ylab=expression(f[X](lambda)),ylim=range(0,max(sigsq))) for(j in 1:k){ lines(rep(lam[which.min((lam+omega[j])^2)],2),c(0,sigsq[j])); points(lam[which.min((lam+omega[j])^2)],sigsq[j],pch=19) lines(rep(lam[which.min((lam-omega[j])^2)],2),c(0,sigsq[j])); points(lam[which.min((lam-omega[j])^2)],sigsq[j],pch=19) } par(mar=c(4,4,2,2)) plot(hvec,gammah,type='p',pch=19,xlab='h',ylab='ACVF',ylim=range(-30,30)) for(j in 1:length(hvec)){lines(c(hvec[j],hvec[j]),c(0,gammah[j]))} @ <>= Uvec<-rnorm(k)*sqrt(sigsq) Vvec<-rnorm(k)*sqrt(sigsq) tvec<-seq(0,500,by=1) Xvec<-tvec*0 for(j in 1:k){ Xvec<-Xvec+Uvec[j]*cos(omega[j]*tvec)+Vvec[j]*sin(omega[j]*tvec) } par(mar=c(4,4,2,2)) plot(tvec,Xvec,type='l',xlab='t',ylab=expression(X[t])) @ Finally suppose $k=50$, with \begin{itemize} \item $\omega_j = j \pi/k$, \item $\sigma_j^2 = 1/j^2$ \end{itemize} for $j=1,\ldots,k$. <>= k<-50 lam<-seq(-pi,pi,by=0.01) omega<-c(1:k)*pi/k sigsq<-1/c(1:k)^2 hvec<-c(0:25) gammah<-hvec*0 for(j in 1:k){ gammah<-gammah+sigsq[j]*cos(omega[j]*hvec) } acvf<-cbind(hvec,gammah) rownames(acvf)<-paste('h=',hvec,sep='') colnames(acvf)<-c('h','ACVF') acvf[1:10,2] @ <>= Flam<-lam*0 par(mar=c(4,4,2,2)) cvec<-c('red','green','blue','black') plot(lam, Flam,type='l',xlab=expression(lambda),ylab=expression(F[X](lambda)),ylim=range(0,sum(sigsq))) for(j in 1:k){ Flamj<-sigsq[j]*((lam >= -omega[j])*0.5+(lam >= omega[j])*0.5) Flam<-Flam+Flamj Flamj[which.min((lam+omega[j])^2)]<-NA Flamj[which.min((lam-omega[j])^2)]<-NA Flam[which.min((lam+omega[j])^2)]<-NA Flam[which.min((lam-omega[j])^2)]<-NA } lines(lam,Flam) for(j in 1:k){ points(lam[which.min((lam+omega[j])^2)+1],Flam[which.min((lam+omega[j])^2)+1],pch=19,cex=0.5) points(lam[which.min((lam-omega[j])^2)-1],Flam[which.min((lam-omega[j])^2)+1],pch=19,cex=0.5) } plot(lam,lam*0,type='n',xlab=expression(lambda),ylab=expression(f[X](lambda)),ylim=range(0,max(sigsq))) for(j in 1:k){ lines(rep(lam[which.min((lam+omega[j])^2)],2),c(0,sigsq[j])); points(lam[which.min((lam+omega[j])^2)],sigsq[j],pch=19) lines(rep(lam[which.min((lam-omega[j])^2)],2),c(0,sigsq[j])); points(lam[which.min((lam-omega[j])^2)],sigsq[j],pch=19) } par(mar=c(4,4,2,2)) plot(hvec,gammah,type='p',pch=19,xlab='h',ylab='ACVF',ylim=range(-2,2)) for(j in 1:length(hvec)){lines(c(hvec[j],hvec[j]),c(0,gammah[j]))} Uvec<-rnorm(k)*sqrt(sigsq) Vvec<-rnorm(k)*sqrt(sigsq) tvec<-seq(0,500,by=1) Xvec<-tvec*0 for(j in 1:k){ Xvec<-Xvec+Uvec[j]*cos(omega[j]*tvec)+Vvec[j]*sin(omega[j]*tvec) } par(mar=c(4,4,2,2)) plot(tvec,Xvec,type='l',xlab='t',ylab=expression(X[t])) @ \pagebreak For these processes, the variance of $\{X_t\}$ can be decomposed into the variances contributed by the harmonic components at different frequencies $\omega_1,\ldots,\omega_k$. In the limit as $k \longrightarrow \infty$ \[ X_t \stackrel{\text{m.s.}}{=} \int_{(-\pi,\pi]} e^{\imagi t \lambda} dZ(\lambda) \] where $\{Z(\lambda): -\pi < \lambda \leq \pi\}$ is a non-decreasing stochastic process. In general, $Z$ is a \emph{complex-valued} process which has uncorrelated increments, taking the form \[ X_t = X_{1,t} + \imagi X_{2,t} \] say, where $\E[X_t] = \E[X_{1,t}] + \imagi \E[X_{2,t}]$ and \[ \E[X_t X_{t+h}] = \E[(X_{1,t} + \imagi X_{2,t})(X_{1,t+h} + \imagi X_{2,t+h})]= \E[X_{1,t} X_{1,t+h} - X_{2,t} X_{2,t+h}] + \imagi \E[X_{1,t} X_{2,t+h} + X_{2,t} X_{1,t+h }] \] etc. We have \[ X_t = \int_{(-\pi,\pi]} e^{\imagi t \lambda} dZ(\lambda) = \int_{(-\pi,\pi]} \left[ \cos (t \lambda) + \imagi \sin (t \lambda) \right] \: dZ(\lambda) \] where $Z(\lambda)$ is a random process on $(-\pi,\pi]$. We need only consider $Z(\lambda)$ for $0 \leq \lambda \leq \pi$. For finite $k$, if \[ dZ (\lambda) = \left\{ \begin{array}{lc} \dfrac{U_j + \imagi V_j}{2} & \lambda = -\omega_j, \: j=1,\ldots,k\\[12pt] \dfrac{U_j - \imagi V_j}{2} & \lambda = \omega_j, \: j=1,\ldots,k\\[12pt] 0 & \text{otherwise} \end{array} \right. \] Here \[ \E[dZ (\lambda) \overline{dZ (\lambda)}] = \E[(U_j^2+V_j^2)/4] = \frac{\sigma_j^2}{2} \] for $\lambda = \pm \omega_j$, $j=1,\ldots,k$ and zero otherwise. We have that \begin{align*} X_t & = \int_{(-\pi,\pi]} e^{\imagi t \lambda} dZ(\lambda) \\[6pt] & = \dfrac{1}{2}\sum_{j=1}^k (U_j + \imagi V_j) e^{-\imagi \omega_j t} + \dfrac{1}{2} \sum_{j=1}^k (U_j - \imagi V_j) e^{\imagi \omega_j t} \\[6pt] & = \sum_{j=1}^k \left[ U_j \cos(\omega_j t) + V_j \sin(\omega_j t) \right] \end{align*} where, by the previous calculation, yields \[ \gamma_X(h) = \sum_{j=1}^k \sigma_j^2 \cos(\omega_j h) \] Recall the earlier inversion formula \[ \gamma_X(h) = \int_{-\pi}^\pi \cos(h \lambda) f_X(\lambda) \: d \lambda = 2 \int_{0}^\pi \cos(h \lambda) f_X(\lambda) \: d \lambda. \] This informs us that for a target $\gamma_X(h)$, we should choose the $\sigma_j^2$ quantities via the spectral density function. We can consider the limit of this process as $k \lra \infty$; we require \[ \sum_{j=1}^\infty \sigma_j^2 < \infty \] which implies (if the $\sigma_j^2$ values are sorted into descending order) that $\sigma_j^2 \lra 0$ as $j \lra \infty$. \bigskip \textbf{Example: } $\AR(1)$ process. \medskip If $X_t \sim \AR(1)$ with parameter $\phi$, $-1 < \phi < 1$, we have that \[ \gamma_X(h) = \frac{\sigma_Z^2}{1-\phi^2} \phi^h \qquad h = 0,1,2\ldots \] and \[ f_X(\lambda) = \frac{\sigma_Z^2}{2 \pi} \frac{1}{|1-\phi e^{-\imagi \lambda}|^2} \qquad -\pi < \lambda < \pi \] This suggests that if we take $\omega_j = j \pi /k, j=0,1,\ldots, k$ as a regular grid with spacing $\pi/k$, we should set \[ \sigma_j^2 = \frac{\sigma_Z^2}{k} \frac{1}{|1-\phi e^{-\imagi \omega_j}|^2} \] as then, if $k$ is large enough \[ \frac{\sigma_0^2}{2} + \sum_{j=1}^{k-1} \sigma_j^2 \cos(h \omega_j) + \frac{\sigma_k^2}{2} \cos(\pi h) \approx \int_{-\pi}^\pi \cos(h \lambda) f_X(\lambda) \ d \lambda \] using a trapezium rule approximation. <>= k<-1000 phi<-0.8 sigZsq<-4 sigZsq/(1-phi^2) lam<-seq(-pi,pi,by=0.01) omega<-c(0:k)*pi/k sigsq<-(sigZsq/k)*(1/Mod(1-phi*exp(-1i*omega))^2) sum(sigsq[2:k])+0.5*(sigsq[1]+sigsq[k+1]) hvec<-c(0:25) gammah<-0.5*sigsq[1]*cos(omega[1]*hvec) for(j in 2:k){ gammah<-gammah+sigsq[j]*cos(omega[j]*hvec) } gammah<-gammah+0.5*sigsq[k+1]*cos(omega[k+1]*hvec) acvf<-cbind(hvec,gammah) rownames(acvf)<-paste('h=',hvec,sep='') colnames(acvf)<-c('h','ACVF') acvf[1:10,2] acvf[1:10,2]/acvf[1,2] #ACF @ <>= Flam<-lam*0 par(mar=c(4,4,2,2)) cvec<-c('red','green','blue','black') plot(lam, Flam,type='n',xlab=expression(lambda),ylab=expression(F[X](lambda)),ylim=range(0,sum(sigsq))) for(j in 1:(k+1)){ Flamj<-sigsq[j]*((lam >= -omega[j])*0.5+(lam >= omega[j])*0.5) Flam<-Flam+Flamj } lines(lam,Flam) plot(omega,k*sigsq/(2*pi),type='l',xlab=expression(lambda),ylab=expression(f[X](lambda)),ylim=range(0,20)) acvf.vals<-sigZsq/(1-phi^2)*phi^hvec lamvec<-seq(0,pi,by=0.001) sdf.vec<-lamvec*0 for(j in 1:length(lamvec)){ sdf.vec[j]<-sum(exp(-1i*lamvec[j]*hvec[-1])*acvf.vals[2:26])+acvf.vals[1]+ sum(exp(1i*lamvec[j]*hvec[-1])*acvf.vals[2:26]) } sdf.vec<-sdf.vec/(2*pi) #lines(lamvec,sdf.vec,col='red') par(mar=c(4,4,2,2)) plot(hvec,gammah/gammah[1],type='p',pch=19,xlab='h',ylab='ACF',ylim=range(-1,1)) for(j in 1:length(hvec)){lines(c(hvec[j],hvec[j]),c(0,gammah[j]/gammah[1]))} Uvec<-rnorm(k+1)*sqrt(sigsq) Vvec<-rnorm(k+1)*sqrt(sigsq) tvec<-seq(0,500,by=1) Xvec<-tvec*0 for(j in 1:(k+1)){ Xvec<-Xvec+Uvec[j]*cos(omega[j]*tvec)+Vvec[j]*sin(omega[j]*tvec) } par(mar=c(4,4,2,2)) plot(tvec,Xvec,type='l',xlab='t',ylab=expression(X[t]),ylim=range(-10,10)) @ \end{document}