课本案例

01章

#安装
install.packages("forecast")

#调用
library(forecast)

#直接录入
#1、行输入
price=c(101,82,66,35,31,7)
  #调用c函数，以行输入的方式将6个时序数据依次赋值给了price这个变量。
  #price以向量的方式存储了这6个数据。
price=ts(price,start=c(2005,1),frequency=12)
  #指定price为时间序列变量，而且序列名仍然取为price。
  #start选项指定序列的起始读入时间。本例指定从2005年1月开始读入。
  #frequency选项指定序列每年读入的数据频率。本例指定读入频率为每年12个，即该序列为月度数据。
  #半年度数据：frequency=2，季度数据：frequency=4，周数据：frequency=52，日数据：frequency=365
price
  #查看时间序列price的情况。运行结果：行标识为年份，列标识为月份，对应的6个月度的时间序列值。

#2、列输入
price=scan()
  #调用scan函数，以列输入的方式将6个时序数据以列的方式读入（手工输入/复制粘贴）
price=ts(price,start=c(2005,1),frequency=12)
price

#外部数据文件转换
x=read.table("E:/data/file1.csv",sep=",",header=T)
  #sep指定制表分隔符。
  #txt格式：sep="\t"：文件以制表符分隔
  #csv格式：sep=","：文件以逗号分隔
  #header指定第一行是否包含变量名。
  #herder=TURE（简写T）：要转换的文本文件第一行包含变量名
  #header=FALSE（简写F）：赶文件第一行不包含文件名
  #x$year或x$yield：“数据框$变量名”：单独考察这两个变量的特征
x
x$Year

#变换
y=log(x$yield)
y
  #时间序列分析得到观察值序列，分析观察值序列的某个函数变换，比如对数变换序列

#子集
z=subset(x,Year>1925,select=yield)
z

#缺失值插值
a=c(1:7)
  #将1~7赋值给变量a
a[4]<NA
  #定义序列a第四个序列值为缺失值
a

#线性插值
y1=na.approx(a)
y1

#样条插值
y2=na.spline(a)
y2
  #线性插值和样条插值都可以由zoo程序包提供。

02章

#时序图
#2.1 默认格式输出
yield=c(15.2,16.9,15.3,14.9,15.7,15.1,16.7)
  #用行输入的方式将7个序列值赋值给向量yield
yield=ts(yield,start=1884)
  #指定yield为时序变量，观察值起始时间为1884年，数据频率为年度数据
plot(yield)
  #绘制yield的时序图，按R语言默认格式输出

#自定义图形参数
#点线结构参数
#type="p"   点     type="o"   线穿过点
#type="l"   线     type="h"    悬垂线
#type="b" 点连线   type="s"    阶梯线

#散点图
plot(yield,type="p")

#点线图
plot(yield,type="o")

#符号参数
plot(yield,type="o",pch=17)

#连线类型参数
#lty=1    实线     lty=4    点+短虚线
#lty=2    虚线     lty=5    长虚线
#lty=3    点线     lty=6    点+长虚线
plot(yield,lty=2)

#线的宽度参数
#lwd=1       默认宽度
#lwd=k     默认宽度的k倍
#lwd=-k   默认宽度的1/k倍
plot(yield,lwd=2)

#颜色参数
#col=1     col="black"     黑色
#col=2     col="red"       红色
#col=3     col="green"     绿色
#col=4     col="blue"      蓝色

#添加文本
plot(yield,main="1884-1890 年英格兰和威尔士地区小麦平均亩产量",xlab="年份",ylab="亩产量")

#指定坐标轴范围
#指定输出横轴范围
plot(yield,xlim=c(1886,1890))

#指定输出纵轴范围
plot(yield,ylim=c(15,16))

#添加参照线
#添加一条垂线
plot(yield)
abline(v=1887,lty=2)

#添加多条垂直参照线
plot(yield)
abline(v=c(1885,1889),lty=2)

#添加水平线
plot(yield)
abline(h=c(15.5,16.5),lty=2)

#绘制序列自相关图
# acf(x,lag=)
#-x:变量名
#-lag:延迟阶数，若用户不特殊指定的话，系统会根据序列长度自动指定延迟阶数

acf(yield)

#2.1 时序图检验
sha=read.table("E:/data/file4.csv",sep=",",header=T)
output=ts(sha$output,start=1964)
plot(output)
#有明显递增趋势→不满足均值、方差为常数→不是平稳时间序列

#2.2 
a=read.table("E:/data/file5.csv",sep=",",header=T)
milk=ts(a$milk,start=c(1962,1),frequency=12)
plot(milk)
#有明显递增趋势及周期性→不是平稳时间序列

#2.3
b=read.table("E:/data/file6.csv",sep=",",header=T)
temp=ts(b$temp,start=1949)
plot(temp)
#最高温度在37度上下波动→是平稳序列

#自相关图检验
#2.1
acf(output,lag=25)
#平面二维悬垂线
#不是平稳时间序列

#2.2
acf(milk)
#不是平稳时间序列

#2.3
acf(temp)
#是平稳序列

#2.4
# rnorm(n=,mean=,sd=)
#n：随机数个数
#mean：均值，缺省值默认为0；
#sd：标准差，缺省值默认为1；
#rnorm函数也可简写为rnorm（n，均值，标准差）
#如果要产生n个服从标准正态分布的随机数，可以简写为rnorm(n)

#标准正态白噪声序列时序图
white_noise=rnorm(1000)
white_noise=ts(white_noise)
plot(white_noise)

#白噪声序列样本自相关图
acf(white_noise)

# Box.test函数
# Box.test(x,type=,lag=)
# X：检验统计量类型
# （1）type="Box-Pierce"，输出白噪声检验的Q统计量，该统计量为系统默认输出结果。
# （2）type="Ljung-Box"，输出白噪声检验的LB统计量。
# -lag：延迟阶数。lag=n 表示输出滞后n阶的白噪声检验统计量，忽略该选项时，默认输出滞后1阶的检验统计量结果。

Box.test(white_noise,lag=6)

Box.test(white_noise,lag=12)

# for函数
# for(x,in,n1:n2) state
# -x：循环变量名
# -n1:n2：给出的循环取值区间
# -state：需要循环执行的命令

for(i in 1:2) print(Box.test(white_noise,lag=6*i))

for(i in 1:2) print(Box.test(temp,lag=6*i))
# 输出有误

for(i in 1:2) print(Box.test(temp,lag=6*i))

03章

04章

05章

#5.1 序列蕴涵显著的线性趋势，1阶差分就可以实现趋势平稳
a=read.table("E:/data/file4.csv",",",header=T)
x=ts(a$output,start=1964)
x.dif=diff(x) 
#diff(x,lag=,differences=)
#x:变量名
#lag:差分的步长，不特意指定，系统默认lag=1
#differences:差分次数，不特意指定，系统默认differences=1
#1阶差分：diff(x)
#2阶差分：diff(x,1,2)
#k阶差分：diff(x,d,k)
#d阶差分：diff(x,d,1),或者简写为diff(x,d)
#一阶差分后再进行d步差分：diff(diff(x),d)
plot(x.dif) #绘制1阶差分序列时序图

#5.2 序列蕴涵曲线趋势，通常低阶（2阶或3阶）差分就可以提取出曲线趋势的影响
b=read.table("E:/data/file15.csv",",",header=T)
x=ts(b$vehicle,start=1950)
plot(x)

#1阶差分，并绘制差分后序列的时序图
x.dif=diff(x)
plot(x.dif)

#2阶差分，并绘制差分后序列时序图
x.dif2=diff(x,1,2)
plot(x.dif2)

#5.3 蕴涵固定周期的序列
c=read.table("E:/data/file16.csv",",",header=T)
x=ts(c$export,start=c(2001,1),frequency=12)
plot(x)

#1阶差分，并绘制差分后序列时序图
x.dif=diff(x)
plot(x.dif)

#1阶差分加12步差分，并绘制差分后时序图
x.dif1_12=diff(diff(x),12)
plot(x.dif1_12)

#5.5 
x=arima.sim(n=1000,list(order=c(0,1,0)),sd=10)
plot(x)

#5.6 
d=read.table("E:/data/file17.csv",",",header=T)
x=ts(d$index,start=1952)
plot(x)

#1阶差分，并绘制差分后序列时序图
x.dif=diff(x)
plot(x.dif)

#绘制差分后序列自相关图和偏自相关图
acf(x.dif) #自相关图
pacf(x.dif) #偏自相关图

#拟合ARIMA(0,1,1)
x.fit=arima(x,order=c(0,1,1))
x.fit

#残差白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

#5.6
library(zoo)
library(forecast)
d=read.table("E:/data/file17.csv",",",header=T)
x=ts(d$index,start=1952)
x.fit=arima(x,order=c(0,1,1))
x.fore=forecast(x.fit,h=10)
x.fore
plot(x.fore)

#arima(x,order=,include.mean=,method=,transform.pars=,fixed=)
#x:要进行模型拟合的序列号
#order:指定模型阶数。order=c(p,d,q),p为自回归阶数，d为差分阶数，q为移动平均阶数
#include.mean:transform.pars=T，需要拟合常数项；transform.pars=F，不需要拟合常数项
#method:指定参数估计方法
#transform.pars:指定参数估计是否由系统自动完成
#transform.pars=T:系统默认设置是系统根据order选项设置的模型阶数自动完成参数估计
#transform.pars=F:需要拟合疏系数模型，不能让系统根据模型的最高阶数自动完成所有参数的估计，我们需要进行人为干扰
#fixed:对疏系数模型指定疏系数的位置

#5.8
e=read.table("E:/data/file18.csv",",",header=T)
x=ts(e$fertility,start=1917)
plot(x)

#1阶差分，并绘制差分后序列时序图
x.dif=diff(x)
plot(x.dif)

#绘制差分后序列自相关图和偏自相关图
acf(x.dif) #自相关图
pacf(x.dif) #偏自相关图

#拟合疏系数模型ARIMA((1,4),1,0)
x.fit=arima(x,order=c(4,1,0),transform.pars=F,fixed=c(NA,0,0,NA))
x.fit

#残差序列白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*1))

#做5期预测，并绘制预测图
x.fore=forecast(x.fit,h=5)
x.fore
plot(x.fore)

#arima(x,order=,include.mean=,method=,transform.pars=,fixed=,seasonal=)
#x:要进行模型拟合的序列名
#order:指定模型阶数
#include.mean:指定是否需要拟合常数项
#method:指定参数估计方法
#transform.pars:指定是否需要人为干扰参数估计
#fixed:对疏系数模型指定疏系数的位置
#seasonal:指定季节模型的阶数与季节周期。该选项的命令格式为：
#seasonal=list(order=c(P,D,Q),period=π)
#加法模型：P=0,Q=0
#乘法模型：P,Q不全为零

#5.9
f=read.table("E:/data/file19.csv",",",header=T)
x=ts(f$unemployment_rate,start=c(1962,1),frequency=4)
plot(x)

#1阶4步差分，并绘制差分后序列时序图
x.dif=diff(diff(x),4)
plot(x.dif)

#差分后序列自相关图和偏自相关图
acf(x.dif) #自相关图
pacf(x.dif) #偏自相关图

#拟合加法季节模型ARIMA((1,4),(1,4),0)
x.fit=arima(x,order=c(4,1,0),seasonal=list(order=c(0,1,0),period=4),transform.par=F,fixed=c(NA,0,0,NA))
x.fit

#残差序列白噪声检验
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

#做3年期预测，并绘制预测图
x.fore=forecast(x.fit,h=12)
x.fore
plot(x.fore)

#5.10
g=read.table("E:/data/file20.csv",sep=",",header=T)
x=ts(g$unemployment_rate,start=c(1948,1),frequency=12)
plot(x)

#作1阶12步差分，并绘制差分后序列时序图
x.dif=diff(diff(x),12)
plot(x.dif)

#绘制差分后序列自相关图和偏自相关图
acf(x.dif) #自相关图
pacf(x.dif) #偏自相关图

#拟合ARIMA(1,(1,12),1)模型
x.fit=arima(x,order=c(1,1,1),seasonal=list(order=c(0,1,0),period=12))
for(i in 1:2) print(Box.test(x.fit$residual,lag=6*i))

#拟合ARIMA(1,1,1)×ARIMA(0,1,1)12模型
x.fit=arima(x,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12))
x.fit

#残差序列白噪声检验
for(i in 1:2) print(Box.test(x.fit$ residual,lag=6*i))

#5.6（续）
#拟合关于时间t的线性回归模型
d=read.table("E:/data/file17.csv",",",header=T)
x=ts(d$index,start=1952)
t=c(1:37)
x.fit1=lm(x~t)
summary(x.fit1)

#拟合关于延迟变量的自回归模型
xlag=x[2:37]
x2=x[1:36]
x.fit2=lm(x2~xlag)
summary(x.fit2)

#两个趋势拟合模型的拟合效果图
fit1=ts(x.fit1$fitted.value,start=1952)
fit2=ts(x.fit2$fitted.value,start=1952)
plot(x,type="p",pch=8)
lines(fit1,col=2)
lines(fit2,col=4)

#Durbin-Watson检验
#检验第一个确定性趋势模型
library(lmtest)
dwtest(x.fit1)

#Durbin检验
#检验第二个确定性趋势模型
dwtest(x.fit2,order.by =xlag)

#残差自相关模型拟合
acf(x.fit1$residual) #残差序列自相关图
pacf(x.fit1$residual) #残差序列偏自相关图

#拟合AR(2)模型
r.fit=arima(x.fit1$residual,order=c(2,0,0),include.mean=F)
r.fit

#残差自相关模型的显著性检验
for(i in 1:2) print(Box.test(r.fit$residual,lag=6*i))

#5.11
h=read.table("E:/data/file21.csv",sep=",",header=T)
x=ts(h$yield_rate,start=c(1963,4),frequency=12)
plot(x)

#作1阶差分，并绘制差分后残差序列时序图
x.dif=diff(x)
plot(x.dif)

#绘制1阶差分后残差平方图
plot(x.dif^2)

#作对数变换，并绘制对数变换后时序图
lnx=log(x)
plot(lnx)

#作1阶差分，并绘制差分后序列时序图
dif.lnx=diff(lnx)
plot(dif.lnx)

#残差白噪声检验
for(i in 1:2) print(Box.test(dif.lnx,lag=6*i))

#ARCH模型 集群效应
#5.12
#时序图
k=read.table("E:/data/file22.csv",sep=",",header=T)
x=ts(k$returns,start=c(1926,1),frequency=12)
plot(x)

#序列平方图
plot(x^2)

#garch(x,order=)
#x:序列名
#order:拟合模型阶数
#拟合ARCH(q)模型为 order=c(0,q)
#拟合GARCH(p,q) 模型为 order=c(p,q)

#LM检验(有误)
for(i in 1:5) print(ArchTest(x,lag=i))

#Portmanteau Q检验
for (i in 1:5) print(Box.test(x^2,lag=i))

#拟合ARCH(3)模型
x.fit=garch(x,order=c(0,3))
summary(x.fit)

#绘制条件异方差模型拟合的95%置信区间
x.pred=predict(x.fit)
plot(x.pred)

#条件异方差置信区间和方差齐性置信区间比较图示
plot(x)
lines(x.pred[,1],col=2)
lines(x.pred[,2],col=2)
abline(h=1.96*sd(x),col=4,lty=2)
abline(h=-1.96*sd(x),col=4,lty=2)

#5.13
w=read.table("E:/data/file23.csv",sep=",",header=T)
x=ts(w$exchange_rates,start=c(1979,12,31),frequency=365)
plot(x)

#P185











#习题1 P193
d=read.table("E:/datali/习题5.1数据.csv",sep=",",header=T)
summary(d)
x=ts(d)
plot(x)
Box.test(x)
x.dif=diff(x);plot(x.dif)
acf(x.dif)
pacf(x.dif)
x.fit=arima(x,order=c(1,1,1));x.fit
library(forecast);librart(zoo)
x.fit2=auto.arima(x);x.fit2
x.fore1=forecast(x.fit);x.fore1
x.fore2=forecast(x.fit);x/fore2
plot(x.fore1,J=283.4);plot(x.fore2)
line(x.fore1)
lines
plot(x.fore1);abline(h=283.4,col=2)

#习题2
#library(forecast);library(zoo)
e=read.table("E:/datali/习题5.2数据.csv",sep=",",header=T)
View(e)
y=ts(e$货运量,start=1949)
plot(y)
y.dif=diff(y);plot(y.dif)
acf(y.dif)
pacf(y.dif)
y.fit1=arima(y,order=c(1,1,1));y.fit1
y.fit2=auto.arima(y);y.fit2
y.dif2=diff(y.dif);plot(y.dif2)
acf(y.dif2)
pacf(y.dif2)

06章

附