Question 4.1.

Part A.

Given the AR2 model \[X_{n}=1.5X_{n-1}-0.8X_{n-2}+\epsilon_{n},\] we have \[\epsilon_{n}=X_{n}-1.5X_{n-1}+0.8X_{n-2}.\] The ACF for \(\left\{\epsilon_{n}\right\}\) is therefore \[\begin{eqnarray} \gamma_{\epsilon}\left(h\right)&=&{\mathrm{Cov}}\left(\epsilon_{n+h},\epsilon_{n}\right)\\ &=&{\mathrm{Cov}}\left(X_{n+h}-1.5X_{n+h-1}+0.8X_{n+h-2},X_{n}-1.5X_{n-1}+0.8X_{n-2}\right)\\ &=&3.89\gamma_{X}\left(h\right)-2.7\left[\gamma_{X}\left(h+1\right)-\gamma_{X}\left(h-1\right)\right]+0.8\left[\gamma_{X}\left(h+2\right)-\gamma_{X}\left(h-2\right)\right]. \end{eqnarray}\]

Since \(\gamma\left(h\right)\) is the inverse Fourier Transformation of the spectral density, we have \[\begin{eqnarray} \gamma\left(h\right)&=&\int_{-\frac{1}{2}}^{\frac{1}{2}}e^{2\pi iwh}f\left(w\right)dw\\ \gamma_{\epsilon}\left(h\right)&=&\int_{-\frac{1}{2}}^{\frac{1}{2}}\left[\left(1+1.5^{2}+0.8^{2}\right)+\left(-1.5-1.5\times0.8\right)\left(e^{2\pi iw}+e^{-2\pi iw}\right)+0.8\left(e^{4\pi iw}+e^{-4\pi iw}\right)\right]e^{2\pi iwh}\lambda_{X}\left(w\right)dw\\ &=&\int_{-\frac{1}{2}}^{\frac{1}{2}}\left[3.89-5.4\cos\left(2\pi w\right)+1.6\cos\left(4\pi w\right)\right]e^{2\pi iwh}\lambda_{X}\left(w\right)dw. \end{eqnarray}\]

The spectrum of the white noise is \[ \lambda_{\epsilon}\left(w\right)=\sum_{h}\gamma_{0}=\sigma^{2}. \] Therefore, by the uniqueness of the Fourier Transformation, we have \[\begin{eqnarray} \lambda_{\epsilon}\left(w\right)&=&\sigma^{2}\\ &=&\left[3.89-5.4\cos\left(2\pi w\right)+1.6\cos\left(4\pi w\right)\right]\lambda_{X}\left(w\right)\\ \Rightarrow\lambda_{X}\left(w\right)&=&\frac{\sigma^{2}}{3.89-5.4\cos\left(2\pi w\right)+1.6\cos\left(4\pi w\right)} \end{eqnarray}\]

The plot of the spectral density and autocovariance are given as follows.

library('TSA')
ARMAspec(model = list(ar = c(1.5, -0.8)))

model_AR = arima.sim(model = list( ar = c(1.5, -0.8)), n=1000)
acf(model_AR, type="covariance")

Part B.

Given the MA2 model

\[ X_{n}=\epsilon_{n-2}+\epsilon_{n-1}+\epsilon_{n} \] with \({\mathrm{Var}}\left(\epsilon_{n}\right)=\sigma^{2}\), we start by writing the ACF \[\begin{eqnarray} \gamma\left(0\right)&=&{\mathrm{Cov}}\left(X_{n},X_{n}\right)=3\sigma^{2}\\ \gamma\left(1\right)&=&\gamma\left(-1\right)={\mathrm{Cov}}\left(X_{n},X_{n-1}\right)=2\sigma^{2}\\ \gamma\left(2\right)&=&\gamma\left(-2\right)={\mathrm{Cov}}\left(X_{n},X_{n-2}\right)=\sigma^{2}\\ \gamma\left(h\right)&=&0\;\forall|h|\geq3. \end{eqnarray}\] So we have \[\begin{eqnarray} \lambda\left(w\right)&=&\sum_{h}\gamma\left(h\right)e^{-2\pi iwh}\\ &=&\sum_{h=-2}^{2}\gamma\left(h\right)e^{-2\pi iwh}\\ &=&\gamma\left(0\right)+\gamma\left(1\right)\cos\left(2\pi w\right)+\gamma\left(2\right)\cos\left(4\pi w\right)\\ &=&\sigma^{2}\left[3+4\cos\left(2\pi w\right)+2\cos\left(4\pi w\right)\right] \end{eqnarray}\] The plot of the spectral density and autocovariance are given as follows:

ARMAspec(model = list(ma = c(1, 1)))

model_MA = arima.sim(model = list( ma = c(1, 1)), n=1000)
acf(model_MA, type="covariance", plot= TRUE)

Part C.

For part A, if we look closely at the spectrum density plot, we can see a peak at frequency around 0.1 (actually slightly less than 0.1), which is the dominant frequency. This indicates that the dominant period is around 10 (precisely larger than 10). If we look at the ACF, we can see there exists an oscillatory behavior characteristic of period = 11. We see that these two are matched with each other.

For part B, there is no appearant peak on the spectrum density plot, and there is no periodic behavior on the ACF plot either. Again, these two correspond to each other.

We can also notice a small oscillation in the spectral density for the MA process, with a dip in the spectral density at frequency 0.33. Thus, for the AR process, the spectral power is concentrated around a peak and the ACF oscillates. For the MA process, the ACF goes quickly to zero and the spectrum oscillates.

Question 4.2.
We first read in the data from the source. We seek to find out the relationship between time and number. We now make a time plot of the data to explore.

mydata = read.table(file="http://ionides.github.io/531w16/hw/sunspots.txt", header = TRUE)
year = mydata$Time
number = mydata$Number
plot(year, number, type= "l")

From the plot, we see that there is a periodic behavior, with regularly spaced peaks. The time interval between these peaks are about 11 years. This is as expected from reading

http://solarscience.msfc.nasa.gov/SunspotCycle.shtml.

Now we take a look at spectrum density.

spectrum(number, main = "Unsmoothed periodogram")

We first use R to see an unparametric method result of the data.

smoothed_r = spectrum(number, spans=c(30,30), main = "Smoothed periodogram")

We now determine the dominant frequency.

smoothed_r$freq[which.max(smoothed_r$spec)]
## [1] 0.007407407

We see that the dominant frequency is 0.007407407, which corresponds to a period of 11.25 years.

Now we use parametric method to estimate the spectral density.

estimated = spectrum(number, method = "ar", main = "Spectrum estimated via AR model picked by AIC")
abline(v=estimated$freq[which.max(estimated$spec)], lty="dotted")

The dominant frequecy is given by

estimated$freq[which.max(estimated$spec)]
## [1] 0.008016032

We see that the dominant frequency is 0.008016032, which corresponds to a period of 10.4 years.

the results from two estimations are very close to each other and they are close to the 11 year period. Therefore we can conclude that both of them are reliable and can be used.