1.

Solutions

a.

#using code given in the book
x1 <- 2*cos(2*pi*1:128*6/100)+3*sin(2*pi*1:128*6/100)
x2 <- 4*cos(2*pi*1:128*10/100)+5*sin(2*pi*1:128*10/100)
x3 <- 6*cos(2*pi*1:128*40/100)+7*sin(2*pi*1:128*40/100)
x <- x1+x2+x3
par(mfrow=c(2,2))
plot.ts(x1,ylim=c(-10,10), main="Series of X_1")
plot.ts(x2,ylim=c(-10,10), main="Series of X_2")
plot.ts(x3,ylim=c(-10,10), main="Series of X_3")
plot.ts(x,ylim=c(-16,16),main="sum")

Although its harder to see for \(\{x_3\}\) it appears that the difference is the fact that these series do not finish the final period, we can see that they begin at a value different then where they end.

b.

From the periodgram below we can see that most important periods occur at \(t=0.4\) and \(t=0.6\). But there appears to be two more spikes both at lower and the high frequencies.

par(mfrow = c(1,1))
P <- abs(2*fft(x)/128)^2; Fr <- 0:127/128
plot(Fr,P,type="o",xlab="Frequency",ylab="Periodogram")

c.

First we plot the time series with the additional white noise.

# using the code given in the book
x1 <- 2*cos(2*pi*1:100*6/100)+3*sin(2*pi*1:100*6/100)
x2 <- 4*cos(2*pi*1:100*10/100)+5*sin(2*pi*1:100*10/100)
x3 <- 6*cos(2*pi*1:100*40/100)+7*sin(2*pi*1:100*40/100)
w <- rnorm(100,0,25)
x <- x1 + x2 + x3 + w
plot.ts(x,main="sum")

# and the periodgram using the code given in the book
P <- abs(2*fft(x)/100)^2; Fr <- 0:99/100
plot(Fr,P,type="o",xlab="Frequency",ylab="Periodogram")

Time Series Plot There appears to be no obvious pattern in the time series. Possibply a trend as we approach a time \(t = 100\)

Periodgram From the periodgram we can see some significant peaks at both the lower and higher frequencies with smaller peaks through out the frequency values.

4.

Solutions

a. First for the white noise series \[E(w_t) = 0 \;\; \text{ and }\;;\ \gamma_w(h) = 1 \;\; for\;\; h = 0\]

For \(x_t\)

\[E(x_t) = E(w_t - \theta w_{t-1}) = E(x_t) - \theta E(w_{t-1}) = 0\] \[\gamma_x(0) = (1 + \theta^2), \;\;and \;\; \gamma_x(1) = \theta, \;\; 0 \;\; otherwise\]

Yes the indeed appear to be stationary as time is not a part of either of the functions above.

b.

To find the power spectrum we simply use the above and apply the formula, so we get

\[f(\omega) = \sum_{-\infty}^\infty \gamma(h)e^{-2\pi i \omega h}\]

\(x_t\) has autocovariance defined only at \(|h| = 0, 1\) so this reduces to

\[f(\omega) = \gamma(0)e^{-2\pi i \omega (0)} + \gamma(1)e^{-2\pi i \omega (1)}\] \[ = \gamma(0) + \gamma(1)e^{-2\pi i \omega (1)}\] \[ = 1 + \theta^2 + \theta e^{-2\pi i \omega}\]

and simplyfing

\[ = 1 + \theta^2 + 2\theta cos(2\pi \omega)\]

8.

Figure 4.31 shows the biyearly smoothed (12-month moving average) num- ber of sunspots from June 1749 to December 1978 with n = 459 points that were taken twice per year; the data are contained in sunspotz. With Exam- ple 4.10 as a guide, perform a periodogram analysis identifying the predom- inant periods and obtaining confidence intervals for the identified periods. Interpret your findings.

Solutions

library(astsa)
# to plot the perdiogram
sp <- spec.pgram(sunspotz,taper=0,log="no")

From the periodgram above we can see that there are two big peaks at about a frquency of .1, further there appears to be a third just over the 0 frequency. Let us investigate this further by finding these values.

# we can get a reverse sorted list of the specs 
sort_specs <- sort(sp$spec, decreasing = TRUE)[c(1,2,4)] # get top 3
#to get the higest spectrum freq peak
#the corresponding frequency
p1 <- sp$freq[sp$spec==sort_specs[1]]; p1
## [1] 0.09166667
#frequency at second highest peak
p2 <- sp$freq[sp$spec==sort_specs[2]]; p2
## [1] 0.1
#third peak
#frequency at third peak
p3 <- sp$freq[sp$spec==sort_specs[3]]; p3
## [1] 0.0125
# when do the cycles occur
cat("cycles are occuring at", 1/p1, 1/p2, 1/p3, "years")
## cycles are occuring at 10.90909 10 80 years

Note here that we can consider cycle of 10 and 10.9 as possibly one, but we continue as seperate instead

Now we will create a the CI for each of these peaks.

make_CI <- function(peak_spec){
  u <- qchisq(0.025,2)
  l <- qchisq(0.975,2)
  c((2*peak_spec)/l,(2*peak_spec)/u)  
}

# for peak 1
make_CI(sort_specs[1])
## [1]    8804.265 1282807.445
# for peak 2 
make_CI(sort_specs[2])
## [1]    8290.672 1207975.406
# for peak 3 
make_CI(sort_specs[3])
## [1]   3428.087 499482.385

Each one of these is far too wide to be of much use for us.

21

Determine the theoretical power spectrum of the series formed by combining the white noise series \(w_t\) to form \[y_t = w_{t-2}+ 4w_{t-1} + 6w_t + 4w_{t+1} + w_{t+2}\] Determine which frequencies are present by plotting the power spectrum.

Solutions

To determine power spectrum we calculate the ACVF to be: \[ \gamma(h) = \begin{cases} 68\sigma_w^2 &|h| = 0\\ 56\sigma_w^2 &|h| = 1\\ 28\sigma_w^2 &|h| = 2\\ 8\sigma_w^2 &|h| = 3\\ \sigma_w^2 &|h| = 4\\ \end{cases} \]

Now we use

\[f(\omega) = \sum_{-\infty}^\infty \gamma(h)e^{-2\pi i \omega h}\]

and pluggin in values to get

\[f(\omega) = \sigma_w^2\left[ 68 + 56e^{-2\pi i \omega} + 28e^{-4\pi i \omega} + 8e^{-6\pi i \omega} + e^{-8\pi i \omega}\right]\]

om <- seq(0,5,0.0001)
power.spec <- 70+112*cos(2*pi*om)+56*cos(4*pi*om)+16*cos(6*pi*om)+2*cos(8*pi*om)
plot(om,power.spec,xlab="frequency",type="l")

It appears to be that all of the frequencies are present.