Notes for STAT 730 HW6  and Lecture                            5/1/2017
=======================

*(II) Preliminary approach: use the roughly quadratic nature of cosine to
approximate the given spectral density by MA(1): X_t = .6*(W_t + (2/3)*W_{t-1})
has spectral density  g(x) = .16*(3.25+3*cos(x)) which over plotted with
f(x) = (1-(x/pi)^2)   gives

> curve(1-(x/pi)^2, 0, pi)
  curve(.36*(13/9+(4/3)*cos(x)),0,pi, add=T, col="blue")

> wvec = rnorm(1001)*.6
  xvec = wvec[-1] +(2/3)*wvec[-1001]

### Now overplot smoothed periodogram with these curves rescaled as densities
###   from (0,pi) to (0,.5) : 

  spec.pgram(xvec, kernel("daniell",4), pad=0, detrend=F, log="n", lty=3)
  curve(.36*(13/9+(4/3)*cos(2*pi*x)),0,.5, add=T, col="blue")
  curve(1-4*x^2, 0, .5, add=T, col="red")
 
## More detailed approach:
## first step: calculate ACF for f(x) = 1-(x/pi)^2 = (4/3)*pi for h=0 and 
##       and  (-1)^h *(-4/pi)/h^2 for h nonzero

### Next find AR process with  approximately this spectral density:

> phi =  acf2AR(c(4*pi/3,rep(c(-1,1),20)*(-4/pi)/(1:40)^2))[40,]
## 40 AR coeffs phi[1:40]

> curve(1/Mod(exp(1i*2*pi*outer(x,0:40))%*%c(1,-phi))^2,0,0.5)
### This gives the desired spectral density, up to constant due to these functions
###  changing the autocovariance to ACF.

## now apply the function ARMAtoMA: 
> thet = ARMAtoMA(ar=phi,ma=NULL, 100)
> curve(Mod(exp(1i*2*pi*outer(x,0:100))%*%c(1,thet))^2,0,0.5,add=T,col="orange")
## essentially perfect overlay
### Note that the variance of this MA process is calculated as
> 1 + sum(thet^2)
[1] 1.224117


> Wt = rnorm(1100)
  Zt = numeric(1000)
  for(i in 1:1000) Zt[i] = sum(c(1,thet)*Wt[(100+i):i])
> var(Wt)
[1] 1.006677
> mean(Zt^2)
[1] 1.217132

> tmp = spec.pgram(Zt, kernel("daniell",10), pad=0, log="n",detrend=F)
### NB actual variance associated with spec density on (0,pi) frequency 
###    interval was 4*pi/3, and spec density on (0,.5) would be 2*pi, so 
###    multiplicative vertical scale on spec density is 1.5*variance, or 
###    here 1.5*1.224 = 1.836

> curve(1.836*(1-4*x^2),0,.5, add=T, col="red")

### THIS IS THE DESIRED OVERLAY PLOT !!!

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

SOME COMMENTS on Box-Pierce-Ljung Test

> Box.test(tmpfit2B$resid, fitdf=11, lag=50)
## n = 403, 
> aux = acf(tmpfit2B$resid, lag.max=50, plot=F)   

> 403*405*sum(aux$acf[2:51]^2/(403-(1:50)))
[1] 32.320    ### not the same as the Box.test default of 30.365, 
           ### but the type="L" or "Ljung-Box" alternative
> round(403*405*aux$acf[2:51]^2/(403-(1:50)), 2)
 [1] 2.47 0.39 0.03 0.00 1.77 0.82 0.03 0.67 0.77 0.85 0.65 1.19 0.01 0.21 0.75
[16] 0.02 2.29 0.04 0.60 0.02 0.32 1.31 0.29 0.00 1.61 3.35 0.48 0.00 0.08 2.22
[31] 1.40 0.04 0.27 1.58 0.14 0.38 0.02 0.27 0.18 0.39 0.20 1.74 0.07 1.46 0.82
[46] 0.01 0.00 0.03 0.01 0.06


>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

COMMENTS on CI's in Spectrogram Plots

>  tmp = spec.pgram(SOI.Rec, k, taper=0)

> plot(tmp, plot.type="coh", ci.lty=2)

  plot(tmp, plot.type="ph", ci.lty=2)

> spec.pgram(SOI,k,taper=0)    ### CI has constant width, plotted in blue
                               ### in upper-right

> spec.pgram(SOI,k,taper=0, log="n")    ### no CI at all