3.7(a)
h = c(0:10)
par(mfcol=c(1,2))
x1 <- (-1.25)^(-h)*(1+.2915*h)
plot(x1,type='h', xlab='h')
abline(h=0)
rho1 <- ARMAacf(ar=c(-1.6,-.64), ma = 0, lag.max = 10)
plot(rho1,type='h', xlab='h')
abline(h=0)
3.7(b)
par(mfcol=c(1,2))
x2 <- 0.8766*(10/9)^(-h)+0.1234*(-2)^(-h)
plot(x2, type='h')
abline(h=0)
rho2 <- ARMAacf(ar=c(.4,.45), ma=0, 10)
plot(rho2,type='h')
3.7(c)
polyroot(c(1, -1.2, 0.85))
## [1] 0.7058824+0.8235294i 0.7058824-0.8235294i
par(mfcol=c(1,2))
x3 <- ((1.08465)^(-h))*1.00241*cos(0.86217*h-0.06939)
plot(x3, type='h')
abline(h=0)
rho3 <- ARMAacf(ar=c(1.2,-0.85), ma=0, 10)
plot(rho3, type='h')
abline(h=0)
3.10
library(astsa)
a <- ts(cmort) #轉換成時間序列資料
acf(a, 48)
(regr = ar.ols(a, order=2, demean=F, intercept=T, na.action=na.pass))
##
## Call:
## ar.ols(x = a, order.max = 2, na.action = na.pass, demean = F, intercept = T)
##
## Coefficients:
## 1 2
## 0.4286 0.4418
##
## Intercept: 11.45 (2.394)
##
## Order selected 2 sigma^2 estimated as 32.32
regr$asy.se.coef #0.04/0.003=10.多 對T值來說夠大
## $x.mean
## [1] 2.393673
##
## $ar
## [1] 0.03979433 0.03976163
resid1 <- regr$resid
plot(resid1)
acf(resid1, na.action=na.pass, main="ACF of residuals")
#resid並非時間序列資料 不能這樣畫
length(cmort)
## [1] 508
regr$ar/regr$asy.se.coef$ar
## , , 1
##
## [,1]
## [1,] 10.77014
## [2,] 11.11090
#sample>500 自由度當作趨近於無限大(常態分配查表)
3(a)
data <- read.csv("C:/Users/WIN/Desktop/Data.csv", header=F, sep=",")
X0 <- data[2,1:314];
X <- as.numeric(X0)
3(b)
par(mfcol=c(1,1))
plot(X, type="l", xlab="week")
3(c)
Y <- log(X)
plot(Y, xlab='week')
qqnorm(Y)#常態分配
hist(Y)
3(d)
fit = lm(Y~time(Y), na.action=NULL)
Z <- diff(Y)
plot(Z, type="o", main="First Difference", xlab='week')
#Z <- resid(fit)
#plot(Z, type="o", main="detrended", xlab='week')
3(e)
par(mfcol=c(2,1))
acf(Z,na.action = na.pass) #sampled ACF
acf(Z,type="partial") #sampled PACF
(regr2 = ar.ols(Z, order=5, demean=F, intercept=T))
##
## Call:
## ar.ols(x = Z, order.max = 5, demean = F, intercept = T)
##
## Coefficients:
## 1 2 3 4 5
## 0.1248 0.2214 0.1235 0.0739 0.0559
##
## Intercept: 0.0009641 (0.008629)
##
## Order selected 5 sigma^2 estimated as 0.02292
regr2$asy.se.coef
## $x.mean
## [1] 0.00862932
##
## $ar
## [1] 0.05619324 0.05350136 0.05487995 0.05400787 0.05131497
3(f)
Z.yw = ar.yw(Z, order=5)
Z.yw$x.mean
## [1] 0.004410287
Z.yw$ar
## [1] 0.07976092 0.21658621 0.04115419 0.13598464 0.08039189
sqrt(diag(Z.yw$asy.var.coef))
## [1] 0.05688829 0.05653995 0.05782775 0.05653995 0.05688829
Z.yw$var.pred
## [1] 0.02853213
3(g)
z <- c(1, -0.07976092, -0.21658621, -0.04115419, -0.13598464, -0.08039189)
polyroot(z)
## [1] 0.336495+1.642245i -1.789820+0.662786i 0.336495-1.642245i
## [4] 1.215130-0.000000i -1.789820-0.662786i
a = polyroot(z)
1/Mod(a)
## [1] 0.5965291 0.5239451 0.5965291 0.8229573 0.5239451
Arg(a)
## [1] 1.368695e+00 2.786941e+00 -1.368695e+00 -3.939967e-14 -2.786941e+00
(2*pi)/Arg(a)
## [1] 4.590641e+00 2.254509e+00 -4.590641e+00 -1.594730e+14 -2.254509e+00
d <- (2*pi)/Arg(a)[1:2]
exp(d)
## [1] 98.557600 9.530615
d
## [1] 4.590641 2.254509
#整理過的資料只是detrend 不需要將算出來的週期再乘上exp 因為資料已經改變 只是不太確定週期單位是多少(週?) 只知道是時間
3(h)
Box.test(z, type="Ljung-Box", lag=Z.yw$order)
##
## Box-Ljung test
##
## data: z
## X-squared = 1.877, df = 5, p-value = 0.8659
3(i)比較並評論
par(mfcol=c(2,1))
ar5_acf <- ARMAacf(Z.yw$ar, ma=0, 28)
acf(Z,28,ylim=range(ar5_acf))
lines(0:28, ar5_acf, col=2, lwd=2) #compare sample ACF with model ACF
ar5_pacf <- ARMAacf(Z.yw$ar, ma=0, 28, pacf=TRUE)
acf(Z,28,type="partial",ylim=range(ar5_pacf))
points(1:28, ar5_pacf, col=2, pch=16) #compare sample PACF with model PACF
