Time series analysis and its applications with R examples（时间序列分析及应用-

高价跪求Time series analysis and its applications with R examples第三版（时间序列分析及应用-R)第四章的答案，或者会做第四章题目的，高价求，跪谢跪谢！可商议价格

Example 4.1
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)
x = x1 + x2 + x3

par(mfrow=c(2,2))
plot.ts(x1, ylim=c(-10,10), main = expression(omega==6/100~~~A^2==13))
plot.ts(x2, ylim=c(-10,10), main = expression(omega==10/100~~~A^2==41))
plot.ts(x3, ylim=c(-10,10), main = expression(omega==40/100~~~A^2==85))
plot.ts(x, ylim=c(-16,16), main="sum")


Example 4.2
P = abs(2*fft(x)/100)^2
Fr = 0:99/100                    
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram")


Example 4.6
par(mfrow=c(3,1))      
arma.spec(log="no", main="White Noise")
arma.spec(ma=.5, log="no", main="Moving Average")  
arma.spec(ar=c(1,-.9), log="no", main="Autoregression")

Example 4.8
x = c(1,2,3,2,1)        
c1 = cos(2*pi*1:5*1/5)
s1 = sin(2*pi*1:5*1/5)
c2 = cos(2*pi*1:5*2/5)
s2 = sin(2*pi*1:5*2/5)
omega1 = cbind(c1, s1)
omega2 = cbind(c2, s2)

anova(lm(x~omega1+omega2))  # ANOVA Table
abs(fft(x))^2/5             # the periodogram (as a check)


Example 4.10
par(mfrow=c(2,1))      
soi.per = mvspec(soi, log="no")                 # using astsa package
# soi.per = spec.pgram(soi, taper=0, log="no")  # using stats package
abline(v=1/4, lty="dotted")
rec.per = mvspec(rec, log="no")
# rec.per = spec.pgram(rec, taper=0, log="no")
abline(v=1/4, lty="dotted")

soi.per$spec[40]  # 0.97223;  soi pgram at freq 1/12 = 40/480
soi.per$spec[10]  # 0.05372;  soi pgram at freq 1/48 = 10/480

# conf intervals -  returned value:
U = qchisq(.025,2)    # 0.05063  
L = qchisq(.975,2)    # 7.37775
2*soi.per$spec[10]/L  # 0.01456
2*soi.per$spec[10]/U  # 2.12220
2*soi.per$spec[40]/L  # 0.26355
2*soi.per$spec[40]/U  # 38.40108

# Repeat lines above using rec in place of soi


Example 4.11
par(mfrow=c(2,1))
k = kernel("daniell",4)
soi.ave = mvspec(soi, k, log="no")                  # using astsa package
# soi.ave = spec.pgram(soi, k, taper=0, log="no")   # using stats package
abline(v=c(.25,1,2,3), lty=2)

soi.ave$bandwidth                 # 0.225; bandwidth, time unit = year
soi.ave$bandwidth/frequency(soi)  # 0.01875; bandwidth, time unit = month

df = soi.ave$df         # df = 16.9875
U = qchisq(.025, df)    # U = 7.555916
L = qchisq(.975, df)    # L = 30.17425
soi.ave$spec[10]        # 0.04952026
soi.ave$spec[40]        # 0.1190800

# intervals
df*soi.ave$spec[10]/L   # 0.02787891
df*soi.ave$spec[10]/U   # 0.1113333
df*soi.ave$spec[40]/L   # 0.06703963
df*soi.ave$spec[40]/U   # 0.2677201

# Repeat above commands with soi replaced by rec, for example:
rec.ave = mvspec(rec, k, log="no")
abline(v=c(.25,1,2,3), lty=2)
# and so on.


Example 4.12
t = seq(0, 1, by=1/200)  # WARNING: using t is bad pRactice because it's reserved- but let's be bad
amps = c(1, .5, .4, .3, .2, .1)
x = matrix(0, 201, 6)
for (j in 1:6) x[,j] = amps[j]*sin(2*pi*t*2*j)
x = ts(cbind(x, rowSums(x)), start=0, deltat=1/200)               
ts.plot(x, lty=c(1:6, 1), lwd=c(rep(1,6), 2), ylab="Sinusoids")
names = c("Fundamental","2nd Harmonic","3rd Harmonic","4th Harmonic","5th Harmonic",
          "6th Harmonic","Formed Signal")
legend("topright", names, lty=c(1:6, 1), lwd=c(rep(1,6), 2))
rm(t)                    # Redemption


Example 4.13
kernel("modified.daniell", c(3,3))          # for a list
plot(kernel("modified.daniell", c(3,3)))    # for a graph

par(mfrow=c(2,1))
k = kernel("modified.daniell", c(3,3))
soi.smo = mvspec(soi, k, log="no")                     # using astsa package
# soi.smo = spec.pgram(soi, k, taper=0, log="no")      # using stats package
  abline(v=1, lty="dotted")   
  abline(v=1/4, lty="dotted")
rec.smo = mvspec(rec, k, log="no")  
# rec.smo = spec.pgram(rec, k, taper=0, log="no")
  abline(v=1, lty="dotted")   
  abline(v=1/4, lty="dotted")

df = soi.smo$df           # df = 17.42618
Lh = 1/sum(k[-k$m:k$m]^2) # Lh = 9.232413
Bw = Lh/480               # Bw = 0.01923419 units = month
soi.smo$bandwidth         # Bw = 0.2308103, units = year (value returned by mvspec)         
soi.smo$bandwidth/12      # Bw = 0.01923419 units = month (value returned by mvspec)  

# An easier way to obtain soi.smo:
soi.smo = mvspec(soi, spans=c(7,7), log="no")                  # using astsa package
# soi.smo = spectrum(soi, spans=c(7,7), taper=0, log="no")     # using stats package


Example 4.14
s0 = mvspec(soi, spans=c(7,7), plot=FALSE)                # using astsa package
s50 = mvspec(soi, spans=c(7,7), taper=.5, plot=FALSE)     
# s0 = spectrum(soi, spans=c(7,7), taper=0, plot=FALSE)   # using stats package
# s50 = spectrum(soi, spans=c(7,7), taper=.5, plot=FALSE)
plot(s0$freq, s0$spec, log="y", type="l", lty=2, ylab="spectrum", xlab="frequency")  # dashed line
lines(s50$freq, s50$spec)   # solid line
title("Tapering")
legend("topright", c("with you", "without you"), lty=1:2)
arrows(1.5,.025,.7,.012, col="purple")
text(1.7,.027, "leakage", col="purple")


Example 4.15
# this part uses Yule-Walker for the estimates, so there's no likelihood involved in this
# I guess you could call it pseudo-AIC

spaic = spec.ar(soi, log="no", ylim=c(0,.3))             # min AIC spec, order = 15
text(frequency(soi)*1/52, .07, substitute(omega==1/52))  # El Nino Cycle
text(frequency(soi)*1/12, .27, substitute(omega==1/12))  # Yearly Cycle
(soi.ar = ar(soi, order.max=30))                         # estimates and AICs
dev.new()
plot(1:30, soi.ar$aic[-1], type="o")                     # plot AICs


# Better comparison of pseudo-ICs
n = length(soi)
AIC = rep(0,30) -> AICc -> BIC
for (k in 1:30){
  fit = ar(soi, order=k, aic=FALSE)
# sigma2 = var(fit$resid, na.rm=TRUE)  # this doesn't seem to work right - use next line
  sigma2 = fit$var.pred                # with this, all IC pick the order 15 model
  BIC[k] = log(sigma2) + (k*log(n)/n)
  AICc[k] = log(sigma2) + ((n+k)/(n-k-2))
  AIC[k] = log(sigma2) + ((n+2*k)/n) }
IC = cbind(AIC, BIC+1)

dev.new()
ts.plot(IC, type="o", xlab="p", ylab="AIC / BIC")
grid()
text(15.2, -1.48, "AIC")
text(15, -1.35, "BIC")


Example 4.18
sr = mvspec(cbind(soi,rec), kernel("daniell",9), plot.type="coh")
# sr = spec.pgram(cbind(soi,rec), kernel("daniell",9), taper=0, plot.type="coh")
sr$df                     # df = 35.8625
f = qf(.999, 2, sr$df-2)  # f = 8.529792
C = f/(18+f)              # C = 0.3188779
abline(h = C)


Example 4.19
par(mfrow=c(3,1))
plot(soi)        # plot data
plot(diff(soi))  # plot first difference
k = kernel("modified.daniell", 6)  # filter weights
plot(soif <- kernapply(soi, k))    # plot 12 month filter
  
dev.new()                          # open new graphics device
mvspec(soif, spans=9, log="no")    # spectral analysis
# spectrum(soif, spans=9, log="no")  
abline(v=12/52, lty="dashed")

dev.new()
w = seq(0, .5, length=500)         # frequency response
FR = abs(1-exp(2i*pi*w))^2
plot(w, FR, type="l", xlab=expression(omega))


Example 4.21
nobs = length(EXP6)  # number of observations
wsize = 256          # window size
overlap = 128        # overlap

ovr = wsize-overlap
nseg = floor(nobs/ovr)-1;            # number of segments
krnl = kernel("daniell", c(1,1))     # kernel
ex.spec = matrix(0, wsize/2, nseg)
for (k in 1:nseg) {
  a = ovr*(k-1)+1
  b = wsize+ovr*(k-1)
  ex.spec[,k] = mvspec(EXP6[a:b], krnl, taper=.5, plot=FALSE)$spec
  }
x = seq(0, 10,  len = nrow(ex.spec)/2)
y = seq(0, ovr*nseg, len = ncol(ex.spec))
z = ex.spec[1:(nrow(ex.spec)/2),]

# below is text version
filled.contour(x,y,log(z),ylab="time",xlab="frequency (Hz)",nlevels=12,col=gray(11:0/11),main="Explosion")

dev.new()  # a nicer version with color
filled.contour(x, y, log(z), ylab="time", xlab="frequency (Hz)", main="Explosion")

dev.new()  # below not shown in text
persp(x,y,z,zlab="Power",xlab="frequency (Hz)",ylab="time",ticktype="detailed",theta=25,d=2,main="Explosion")  


Example 4.22
library(wavethresh)
eq = scale(EQ5)  # standardize the series
ex = scale(EXP6)

eq.dwt = wd(eq, filter.number=8)
ex.dwt = wd(ex, filter.number=8)

# plot the wavelet transforms
par(mfrow = c(1,2))
plot(eq.dwt, main="Earthquake")   
plot(ex.dwt, main="Explosion")  

# total power
TPe = rep(NA,11)  # for the earthquake series
for (i in 0:10){TPe[i+1] = sum(accessD(eq.dwt, level=i)^2)}
TotEq = sum(TPe)  # check with sum(eq^2)
TPx = rep(NA,11)  # for the explosion series
for (i in 0:10){TPx[i+1] = sum(accessD(ex.dwt, level=i)^2)}
TotEx = sum(TPx)  # check with sum(ex^2)

# make a nice table
Power = round(cbind(11:1, 0:10, 100*TPe/TotEq, 100*TPx/TotEx), digits=3)
colnames(Power) = c("Level", "Resolution", "EQ(%)", "EXP(%)")
Power


Example 4.23
library(wavethresh)
eq = scale(EQ5)
par(mfrow=c(3,1))
eq.dwt = wd(eq, filter.number=8)
eq.smo = wr(threshold(eq.dwt, levels=5:10))
ts.plot(eq, main="Earthquake", ylab="Data")
ts.plot(eq.smo, ylab="Signal")
ts.plot(eq-eq.smo, ylab="Resid")


Example 4.24
LagReg(soi, rec, L=15, M=32, threshold=6)
LagReg(rec, soi, L=15, M=32, inverse=TRUE, threshold=.01)


Example 4.25
SigExtract(soi, L=9, M=64, max.freq=.05)


Example 4.26
per = abs(fft(soiltemp-mean(soiltemp))/sqrt(64*36))^2      
per2 = cbind(per[1:32,18:2], per[1:32,1:18])   # this and line below is just rearranging
per3 = rbind(per2[32:2,], per2)                # results to get 0 frequency in the middle

par(mar=c(1,2.5,0,0)+.1)
persp(-31:31/64, -17:17/36, per3, phi=30, theta=30, expand=.6, ticktype="detailed", xlab="cycles/row",
       ylab="cycles/column", zlab="Periodogram Ordinate")