setwd(“C:\Users\TDL\Desktop\观影人数时间序列”)
library(MASS)
library(forecast)
source(“trndseas.R”)
data=read.table(“chicago.txt”,header=F)
head(data)
names(data)=c(“每日平均收入”,“年月日”)
data
年
月
日
=
a
s
.
D
a
t
e
(
d
a
t
a
年月日=as.Date(data
年月日=as.Date(data年月日,"%m/%d/%Y")
data
星
期
=
w
e
e
k
d
a
y
s
(
d
a
t
a
星期=weekdays(data
星期=weekdays(data年月日)
data$n=1:nrow(data)
######数据趋势
y=c(dataKaTeX parse error: Expected 'EOF', got '#' at position 19: …均收入) ts.plot(y)#̲从图中可以看出数据有明显的周期…每日平均收入^lambda-1)/lambda
lm.lam<-lm(Ylam~dataKaTeX parse error: Expected 'EOF', got '#' at position 20: …summary(lm.lam)#̲box-cox后的模型R方明显…coefficients[1]
beta1<-lm.lam
c
o
e
f
f
i
c
i
e
n
t
s
[
2
]
c
u
r
v
e
(
(
1
+
l
a
m
b
d
a
∗
(
b
e
t
a
0
+
b
e
t
a
1
∗
x
)
)
(
1
/
l
a
m
b
d
a
)
,
f
r
o
m
=
m
i
n
(
d
a
t
a
coefficients[2] curve((1+lambda*(beta0+beta1*x))^(1/lambda), from=min(data
coefficients[2]curve((1+lambda∗(beta0+beta1∗x))(1/lambda),from=min(datan), to=max(data
n
)
,
c
o
l
=
"
b
l
u
e
"
,
l
w
d
=
2
,
x
l
a
b
=
"
n
"
,
y
l
a
b
=
"
每
日
平
均
收
入
"
)
p
o
i
n
t
s
(
d
a
t
a
n), col="blue", lwd=2, xlab="n", ylab="每日平均收入") points(data
n),col="blue",lwd=2,xlab="n",ylab="每日平均收入")points(datan,data$每日平均收入, pch=21, cex=1.2, col=“red”, bg=“orange”)
mtext(“Box-Cox Transformations”, outer = TRUE, cex=1.5)
par(op)
###第四张图给出曲线图和对应的散点图
######ARMA模型
#1.画出时序图和自相关图
ts.plot(y)#从图中可以看出数据有明显的周期性、趋势性
acf(y,main=“自相关图”)#利用acf()函数画出序列的自相关图,通过自相关图判断序列是否平稳
#2.对序列进行平稳化处理
y=log(y)
par(mfrow=c(2,1))
plot(y,type=“l”,xlab=“n”,ylab=“每日平均收入”,main=“差分前”)
acf(y,main=“自相关图”,xlab=“滞后阶数”)
ndiffs(y)#结果表明序列需要进行1阶差分
ndata<-diff(y,1)
ndiffs(ndata)#结果表明无需进行
par(mfrow=c(2,1))
plot(ndata,type=“l”,main=“差分后”)
acf(ndata,main=“自相关图”,xlab=“滞后阶数”)#从上图可以看到,1阶差分以后序列变为平稳序列,且自相关图显示自相关系数在滞后1阶后就快速的减为0。进一步表明序列平稳。
#3、ARMA模型的定阶及参数估计
par(mfrow=c(2,1))
acf(ndata)
pacf(ndata)
model1<-arima(y,order=c(2,1,2),method=“ML”)
model1
#4、模型检验
qqnorm(model1
r
e
s
i
d
u
a
l
s
)
q
q
l
i
n
e
(
m
o
d
e
l
1
residuals) qqline(model1
residuals)qqline(model1residuals)
Box.test(model1KaTeX parse error: Expected 'EOF', got '#' at position 29: …e="Ljung-Box") #̲5、预测 p=predict(…pred)