ts_exp4
李馨雨
2022-05-15
（1） 模型识别 主要通过序列的自相关函数、偏自相关函数表现的特征，进行初步的模型识别.
（2） 模型参数估计
a. 估计MA模型参数
b.
估计模型以后要能写出模型的形式(差分方程形式和用滞后算子表示的形式)
（3） 模型的诊断检验 a.
根据模型残差是不是白噪声来判断模型是否为适应性模型
b. 能根据输出结果判断模型是否平稳，是否可逆
c. 若有多个序列是模型的适应性模型，会用合适的方法从这些模型中进行选择，如比较模型的残差方差，AIC，SC等。
#####先提供几个自定义R函数。
#######绘制MA理论谱密度曲线的函数，模拟MA序列的函数， 从自协方差函数利用递推方法估计模型参数的函数：
## 给定MA系数MA的理论谱密度
ma.true.spectrum <- function(a, ngrid=256, sigma=1,
tit="True MA Spectral Density",
plot.it=TRUE){
p <- length(a)
freqs <- seq(from=0, to=pi, length=ngrid)
spec <- numeric(ngrid)
for(ii in seq(ngrid)){
spec[ii] <- 1 + sum(complex(mod=a, arg=freqs[ii]*seq(p)))
}
spec = sigma^2 / (2*pi) * abs(spec)^2
if(plot.it){
plot(freqs, spec, type='l',
main=tit,
xlab="frequency", ylab="spectrum",
axes=FALSE)
axis(2)
axis(1, at=(0:6)/6*pi,
labels=c(0, expression(pi/6),
expression(pi/3), expression(pi/2),
expression(2*pi/3), expression(5*pi/6), expression(pi)))
box()
}
invisible(list(frequencies=freqs, spectrum=spec,
ma.coefficients=a, sigma=sigma))
}
## 模拟 MA(q) 模型
## a =( a_1, \dots, a_q)'是个列向量
## 模型是 X_t = epsilon_t + a_1 epsilon_{t-1} + \dots + a_q epsilon_{t-q}
## \Var(\epsilon_t) = sigma^2
## 我们使用的是filter函数
ma.gen <- function(n, a, sigma=1.0, by.roots=FALSE,
plot.it=FALSE){
if(by.roots){
require(polynom)
cf <- Re(c(poly.calc(a)))
cf <- cf / cf[1]
a <- cf
}
q <- length(a)
n2 <- n+q
eps <- rnorm(n2, 0, sigma)
x2 <- filter(eps, c(1,a), method="convolution", side=1)
x <- x2[(q+1):n2]
x <- ts(x)
attr(x, "model") <- "MA"
attr(x, "b") <- a
attr(x, "sigma") <- sigma
if(plot.it) plot(x)
x
}
## 给定 \gamma_0, \gamma_1, \dots, \gamma_q,
## 解出MA(q) 系数 b_1, \dots, b_q, \sigma^2
##利用李雷(1991)的算法, 课本P93
## 输入: gms -- \gamma_0, \gamma_1, \dots, \gamma_q
ma.solve <- function(gms, k=100){
q <- length(gms)-1
if(q==1){
rho1 <- gms[2] / gms[1]
b <- (1 - sqrt(1 - 4*rho1^2))/(2*rho1)
s2 <- gms[1] / (1 + b^2)
return(list(b=b, s2=s2))
}
A <- matrix(0, nrow=q, ncol=q)
for(j in seq(2,q)){
A[j-1,j] <- 1
}
cc <- numeric(q); cc[1] <- 1
gamma0 <- gms[1]
gammas <- numeric(q+k)
gammas[1:(q+1)] <- gms
gamq <- gms[-1]
Gammak <- matrix(0, nrow=k, ncol=k)
for(ii in seq(k)){
for(jj in seq(k)){
Gammak[ii,jj] <- gammas[abs(ii-jj)+1]
}
}
Omk <- matrix(0, nrow=q, ncol=k)
for(ii in seq(q)){
for(jj in seq(k)){
Omk[ii,jj] <- gammas[ii+jj-1+1]
}
}
PI <- Omk %*% solve(Gammak, t(Omk))
s2 <- gamma0 - c(t(cc) %*% PI %*% cc)
b <- 1/s2 * c(gamq - A %*% PI %*% cc)
return(list(b=b, s2=s2))
}arma.Wold <- function(n, a, b=numeric(0)){
p <- length(a)
q <- length(b)
arev <- rev(a) #逆序
psi <- numeric(n)
psi[1] <- 1
for(j in seq(n-1)){
if(j <= q) bj=b[j]
else bj=0
psis <- psi[max(1, j+1-p):j]
np <- length(psis)
if(np < p) psis <- c(rep(0,p-np), psis)
psi[j+1] <- bj + sum(arev * psis)
}
psi
}
## Calculate theoretical autocovariance function of ARMA model using Wold expansion
arma.gamma.by.Wold <- function(n, a, b=numeric(0), sigma=1){
nn <- n + 100
psi <- arma.Wold(nn, a, b)
#print(psi)
gam <- numeric(n)
for(ii in seq(0, n-1)){
gam[ii+1] <- sum(psi[1:(nn-ii)] * psi[(ii+1):nn])
}
gam <- (sigma^2) * gam
gam
}
arma.gamma <- arma.gamma.by.Wold第一题
模拟系数分别为0.4、-0.4和5/2的MA(1)模型(n=500),分别画出时序图和acf与pacf图. 并画出系数为0.4与-0.4时的理论谱密度和模拟数据的谱密度,比较一下有什么区别.
#估计协方差函数
gamma<-function(h,n,X){ #传入样本个数n，计算gamma(h)_hat
gammah<-0
for(i in 1:(n-h)){
gammah<-gammah+(X[i]-mean(X))*(X[i+h]-mean(X))
}
return(gammah/n)
}library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(fUnitRoots)## 载入需要的程辑包：timeDate## 载入需要的程辑包：timeSeries## 载入需要的程辑包：fBasics#生成三个系数的模拟序列
#set.seed(1234)
X1.T1<-ma.gen(n=500, a=c(0.4), sigma=1.0, by.roots=FALSE,plot.it=FALSE)
X2.T1<-ma.gen(n=500, a=c(-0.4), sigma=1.0,by.roots=FALSE,plot.it=FALSE)
X3.T1<-ma.gen(n=500, a=c(5/2), sigma=1.0, by.roots=FALSE,plot.it=FALSE)
plot(X1.T1,ylab="",main="MA(1)时序图");acf(X1.T1);pacf(X1.T1)
plot(X2.T1,ylab="",main="MA(1)时序图");acf(X2.T1);pacf(X2.T1)
plot(X3.T1,ylab="",main="MA(1)时序图");acf(X3.T1);pacf(X3.T1)
#adf.test(X3.T1)
#系数为0.4与-0.4的理论谱密度与模拟数据谱密度
plot1.spec<-function(a,X){ #输入MA(q)真实的系数，输入序列
n<-length(X)
q=length(a)
lams <- seq(0, pi, length=201)
spec0 <- ma.true.spectrum(a,ngrid=length(lams),sigma=1,plot.it=FALSE)$spectrum
plot(lams, spec0,lwd=1,type="l",lty=1,main="",xlab="frequency", ylab="")
gams<-c(gamma(0,n,X))
for(i in 1:q){ #生成模拟数据估计自协方差
gams<-c(gams,gamma(i,n,X))
}
res<-ma.solve(gams) #b是模型参数，s2是估计方差
spec1 <- ma.true.spectrum(res$b,ngrid=length(lams),sigma=sqrt(res$s2),plot.it=FALSE)$spectrum
lines(lams, spec1, col="green", lwd=2, lty=2)
legend(0,0.2, lwd=c(1,2), lty=c(1,2),col=c("black", "green"),
legend=c("True","估计谱密度"))
}
plot1.spec(a=c(0.4),X1.T1)
title(main="MA(1)序列谱密度图(theta1=0.4)")
plot1.spec(a=c(-0.4),X2.T1)
title(main="MA(1)序列谱密度图(theta1=-0.4)")
从三张自相关图拖尾和偏相关图截尾,以及做单位根检验均可以看出三个序列都是平稳序列.
低阶MA模型的谱都很平缓。理论谱密度和模拟数据谱密度相差不大，二者均有当系数为b1=0.4>0时谱密度峰出现在低频，b1=-0.4<0时谱密度峰出现在高频的特点。
第二题
模拟如下MA（2）模型:
\(X_t=\epsilon_t-2/5\epsilon_{t-1}+4/25\epsilon_{t-2},\epsilon_t
i.i.d~N(0,1)\)
n=200.
同样的画出模拟数据的acf和pacf图以及理论谱密度和模拟数据的谱密度.
set.seed(12) #复现结果
X1.T2<-ma.gen(n=200,a=c(-2/5,4/25),sigma=1.0,by.roots=FALSE,plot.it=FALSE)
plot(X1.T2,ylab="",main="MA(2)时序图");acf(X1.T2);pacf(X1.T2)
plot1.spec(a=c(-2/5,4/25),X1.T2)
title(main="MA(2)序列谱密度图(theta1=-0.4,theta2=0.16)")
第三题
模拟n=500的ARMA序列：\(X_t=0.8X_{t-1}+\epsilon_t-0.6\epsilon_{t-1},\epsilon_t~iid~N(0,1)\).画出时序图，并根据模拟的序列建模. 利用auto.arima.
library(forecast)
library(tseries)
#模拟数据
X.T3 <- arima.sim(model=list(ar=c(0.8),ma=c(-0.6)), n=500)
plot(X.T3,main="ARMA(1,1) series",ylab="")
#平稳性检验
par(mfrow=c(2,1))
acf(X.T3);pacf(X.T3)
adf.test(X.T3)## Warning in adf.test(X.T3): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: X.T3
## Dickey-Fuller = -6.5354, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary#白噪声检验
Box.test(X.T3,type="Ljung-Box",lag=6)##
## Box-Ljung test
##
## data: X.T3
## X-squared = 106.48, df = 6, p-value < 2.2e-16Box.test(X.T3,type="Ljung-Box",lag=12)##
## Box-Ljung test
##
## data: X.T3
## X-squared = 108.88, df = 12, p-value < 2.2e-16从ACF和PACF图中可以看出序列是平稳序列，单位根检验也验证了这一点。
白噪声检验显示序列值彼此之间蕴涵着相关关系，为非白噪声序列。
#模型定阶
fit1 = auto.arima(X.T3,ic="aic")
fit2 = auto.arima(X.T3,ic="bic")
fit3 = arima(X.T3, order = c(1, 0, 1),method="ML") #真实阶数
fit1;fit2;fit3## Series: X.T3
## ARIMA(1,0,1) with zero mean
##
## Coefficients:
## ar1 ma1
## 0.8004 -0.5790
## s.e. 0.0579 0.0769
##
## sigma^2 = 0.8984: log likelihood = -681.79
## AIC=1369.57 AICc=1369.62 BIC=1382.22## Series: X.T3
## ARIMA(1,0,1) with zero mean
##
## Coefficients:
## ar1 ma1
## 0.8004 -0.5790
## s.e. 0.0579 0.0769
##
## sigma^2 = 0.8984: log likelihood = -681.79
## AIC=1369.57 AICc=1369.62 BIC=1382.22##
## Call:
## arima(x = X.T3, order = c(1, 0, 1), method = "ML")
##
## Coefficients:
## ar1 ma1 intercept
## 0.7917 -0.5726 -0.1011
## s.e. 0.0599 0.0784 0.0863
##
## sigma^2 estimated as 0.8925: log likelihood = -681.11, aic = 1370.23可以看出按照auto.arima函数，AIC准则和BIC准则均给出的是ARMA(1,1)模型，真实阶数模型结合模拟数据给出的系数估计与前两个给出的差别不大.为了进一步验证这一点，对建好的模型做残差诊断。
#残差诊断
par(mfrow=c(1,3))
plot(fit1$residuals, type = "p", pch = "+")
plot(fit2$residuals, type = "p", pch = "+")
plot(fit3$residuals, type = "p", pch = "+")
qqnorm(fit1$residuals)
qqnorm(fit2$residuals)
qqnorm(fit3$residuals)
对三个模型分别做 qq图(考虑正态性)与散点图(考虑随机性)，发现两个模型的残差都非常接近白噪声.
Box.test(fit1$residuals,type="Ljung-Box",lag=log(length(fit1$residuals)))##
## Box-Ljung test
##
## data: fit1$residuals
## X-squared = 1.8864, df = 6.2146, p-value = 0.9396Box.test(fit2$residuals,type="Ljung-Box",lag=log(length(fit2$residuals)))##
## Box-Ljung test
##
## data: fit2$residuals
## X-squared = 1.8864, df = 6.2146, p-value = 0.9396Box.test(fit3$residuals,type="Ljung-Box",lag=log(length(fit3$residuals)))##
## Box-Ljung test
##
## data: fit3$residuals
## X-squared = 1.9612, df = 6.2146, p-value = 0.9336在白噪声检验下，确实是
ARMA(1,1)的残差结果要更好一些,与前面的结果一致， 所以最后我们选择
ARMA(1,1)模型完成数据的建模。
因此拟合的模型为 \(X_t+0.1013=\frac{(1-0.5726B)\epsilon_t}{1-0.7916B},Var(\epsilon_t)=0.8925\).
(每次运行因为模拟序列的随机性会导致系数有微小差异)
第四题
模拟不同参数的ARMA（1,1）过程：$ (,)=(0.4,0.5),(-0.9,-0.5)$ 并画出时序图和acf，pacf. 对比并做比较。利用arima.sim
#模拟数据
X1.T4 <- arima.sim(model=list(ar=c(0.4),ma=c(0.5)), n=500)
X2.T4 <- arima.sim(model=list(ar=c(-0.9),ma=c(-0.5)), n=500)
#平稳性检验
par(mfrow=c(3,1),mar=c(2,2,0,0),mgp=c(3, 1, 0),oma=c(0,0,2,0))
tits <- c("ARMA(1,1) series with phi=0.4,theta=0.5","ARMA(1,1) series with phi=-0.9,theta=-0.5")
plot(X1.T4,main="",ylab="")
acf(X1.T4);pacf(X1.T4)
plot(X2.T4,main="",ylab="")
acf(X2.T4);pacf(X2.T4)
mtext(tits[2],side=3, outer=TRUE)
\((\phi,\theta)=(0.4,0.5)\)这组模型的时序图是连续的，没有趋势特征，ACF减小速度很快，在2阶后开始拖尾，PACF图减小速度很快，在2阶后开始截尾。而\((\phi,\theta)=(-0.9,-0.5)\)这组模型的时序图是间断的（数据骤变次数较多），没有趋势特征，ACF减小速度较慢（比第一个模型显著慢了很多），在25阶后开始拖尾，PACF图减小速度较快，在3阶后开始截尾。
第五题
现有201个连续的生产记录，数据见文件“3.1.txt”（行数据）。
（1）画出时序图，并判断该序列的平稳性和纯随机性；
（2）如果序列平稳且非白噪声，选择适当的模型拟合该序列；
（3）利用拟合模型，预测该序列下一时刻95%的置信区间。
（1）
data<-read.table("3.1.txt",header=FALSE,na.strings = "")
X.T5<- as.matrix(data);X.T5<-na.omit(as.vector(t(X.T5)))
X.T5<-ts(X.T5)
plot(X.T5,main="201个连续生产记录时序图",ylab="")
#平稳性
acf(X.T5);pacf(X.T5)
adf.test(X.T5)##
## Augmented Dickey-Fuller Test
##
## data: X.T5
## Dickey-Fuller = -6.324, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary#随机性
Box.test(X.T5,type="Ljung-Box",lag=log(length(X.T5)))##
## Box-Ljung test
##
## data: X.T5
## X-squared = 25.866, df = 5.3033, p-value = 0.0001264Box.test(X.T5,type="Ljung-Box",lag=6)##
## Box-Ljung test
##
## data: X.T5
## X-squared = 31.083, df = 6, p-value = 2.444e-05Box.test(X.T5,type="Ljung-Box",lag=12)##
## Box-Ljung test
##
## data: X.T5
## X-squared = 33.962, df = 12, p-value = 0.0006838由acf图和pacf图及单位根检验得出序列是平稳的，由白噪声检验得出序列是非白噪声序列。
（2）选择适当的模型拟合该序列
根据自相关图和偏相关图，选择ARMA(3,0)和ARMA(0,1)进行拟合，以及 foreca st 自带的 auto.arima 函数进行建模,并比较.
fit1<-arima(X.T5,order=c(3,0,0),method="ML")
fit2<-arima(X.T5,order=c(0,0,1),method="ML")
fit3<-auto.arima(X.T5,ic="aic")
fit4<-auto.arima(X.T5,ic="bic")
fit1;cat("\n")##
## Call:
## arima(x = X.T5, order = c(3, 0, 0), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 intercept
## -0.4118 -0.2963 -0.1887 84.1273
## s.e. 0.0695 0.0722 0.0701 0.0987
##
## sigma^2 estimated as 6.989: log likelihood = -480.77, aic = 971.53fit2;cat("\n")##
## Call:
## arima(x = X.T5, order = c(0, 0, 1), method = "ML")
##
## Coefficients:
## ma1 intercept
## -0.4775 84.1297
## s.e. 0.0670 0.0991
##
## sigma^2 estimated as 7.163: log likelihood = -483.21, aic = 972.42fit3;cat("\n")## Series: X.T5
## ARIMA(0,0,2) with non-zero mean
##
## Coefficients:
## ma1 ma2 mean
## -0.4106 -0.1070 84.1308
## s.e. 0.0707 0.0652 0.0911
##
## sigma^2 = 7.177: log likelihood = -481.9
## AIC=971.81 AICc=972.01 BIC=985.02fit4## Series: X.T5
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 mean
## -0.4775 84.1297
## s.e. 0.0670 0.0991
##
## sigma^2 = 7.235: log likelihood = -483.21
## AIC=972.42 AICc=972.54 BIC=982.33可以看出选择 AIC准则时，auto.arima给出的是ARMA(0,2),选择BIC准则时，auto.arima给出的是ARMA(0,1),综合各项估计数据,也是ARMA(0,1)会稍微好一些，但差别不大.为了进一步验证这一点，对建好的模型做残差诊断.
#残差诊断
fit3<-arima(X.T5,order=c(0,0,2),method="ML")
par(mfrow=c(1,3))
plot(fit1$residuals, type = "p", pch = "+")
plot(fit2$residuals, type = "p", pch = "+")
plot(fit3$residuals, type = "p", pch = "+")
qqnorm(fit1$residuals)
qqnorm(fit2$residuals)
qqnorm(fit3$residuals)
对三个模型分别做 qq图(考虑正态性)与散点图(考虑随机性)，发现两个模型的残差都非常接近白噪声.
Box.test(fit1$residuals,type="Ljung-Box",lag=log(length(fit1$residuals)))##
## Box-Ljung test
##
## data: fit1$residuals
## X-squared = 0.36135, df = 5.3033, p-value = 0.9976Box.test(fit2$residuals,type="Ljung-Box",lag=log(length(fit2$residuals)))##
## Box-Ljung test
##
## data: fit2$residuals
## X-squared = 5.4411, df = 5.3033, p-value = 0.4022Box.test(fit3$residuals,type="Ljung-Box",lag=log(length(fit3$residuals)))##
## Box-Ljung test
##
## data: fit3$residuals
## X-squared = 2.7033, df = 5.3033, p-value = 0.7797Box.test(fit1$residuals,type="Ljung-Box",lag=12)##
## Box-Ljung test
##
## data: fit1$residuals
## X-squared = 7.9192, df = 12, p-value = 0.7914Box.test(fit2$residuals,type="Ljung-Box",lag=12)##
## Box-Ljung test
##
## data: fit2$residuals
## X-squared = 11.868, df = 12, p-value = 0.4563Box.test(fit3$residuals,type="Ljung-Box",lag=12)##
## Box-Ljung test
##
## data: fit3$residuals
## X-squared = 9.6771, df = 12, p-value = 0.6443在白噪声检验下，确实是 ARMA(3,0)的残差结果要更好一些,所以最后我们选择 ARMA(3,0)模型完成数据的建模.
#参数显著性检验
fit1##
## Call:
## arima(x = X.T5, order = c(3, 0, 0), method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 intercept
## -0.4118 -0.2963 -0.1887 84.1273
## s.e. 0.0695 0.0722 0.0701 0.0987
##
## sigma^2 estimated as 6.989: log likelihood = -480.77, aic = 971.53# ar1 系数的显著性检验
t1 = -0.4118/0.0695
# 利用 P 值函数 pt()计算P值，其中自由度df为时间序列样本数减去未知参数个数，本案例中为 201-4=197
p1 = 2*pt(t1, df=197, lower.tail = T)
# ar2 系数的显著性检验
t2 = -0.2963/0.0722
p2 = 2*pt(t2, df=197, lower.tail = T)
# ar3 系数的显著性检验
t3 = -0.1887/0.0701
p3 = 2*pt(t3, df=197, lower.tail = T)
# 常数的显著性检验
t0 = 84.1273/0.0987
p0 = 2*pt(t0, df=197, lower.tail = F)
p_value = data.frame(p1, p2,p3, p0)
row.names(p_value) = "P 值"
p_value## p1 p2 p3 p0
## P 值 1.371831e-08 5.945963e-05 0.007716623 0检验结果显示，四个系数的P值都小于\(\alpha=0.05\)，拒绝系数等于0的原假设，因此在置信水平
95% 下均显著非零。
因此拟合的模型为\(X_t-84.1273=\frac{\epsilon_t}{1+0.4118B+0.2963B^2+0.1887B^3},Var(\epsilon_t)=6.989\).
（3）利用拟合模型，预测该序列下一时刻95%的置信区间。
x.fore = forecast(fit1, h=1)
x.fore## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 202 85.32682 81.93884 88.7148 80.14535 90.50829所以95%的置信区间为[80.14535,90.50829]
# 预测图
plot(x.fore)