R语言之时间序列图绘制

01 读取数据

02 赋值变量

03 时间序列值绘制

04 读取包

05 单位根检验，又称平稳性检验

06 自相关图和偏相关图

07 白噪声检验或纯随机性检验

08 模型拟合

09 显著性检验

10 预测未来趋势

结语：

我使用的是：R+Rstudio

01 读取数据

读取数据（read.table），数据位置("F:/笔记/Rstudio/A1_7.csv")

首先找到文件

点击上方文件夹获取文件位置

复制文件位置“F:\笔记\Rstudio”后将“\”改为“/”，在后面加上文件名

改变后："F:/笔记/Rstudio/A1_7.csv"

sep=","（空格位置用“，”代替）

header=T（数据第一行是标题，不属于数据）

data（将数据赋值给data）

data<-read.table("F:/笔记/Rstudio/A1_7.csv",sep=",",header=T)   #data数据

02 赋值变量

data$number（取数据data中的名字为number的一列，依据数据具体名称来定，是需要研究的部分）

将提取的一列变量"x"

在这个窗口可以看见

x<-ts(data$number,start=1900)
#ts时间序列   x时间序列未知量

03 时间序列值绘制

绘制时间序列图用plot()

研究变量为：x

main：主标题

xlab：x轴变量名称

ylab：y轴变量名称

xlim：x轴取值范围

ylim：y轴取值范围

plot(x,main = "1900-1998年全球7级以上地震发生次数",xlab = "时间",ylab="7级以上地震发生次数",xlim = c(1900,1998),ylim=c(0,50)
)

绘图如下：

可以通过col=改变线的颜色，可以是1、2、3、4代表黑、红、绿、蓝，或是直接打出颜色名称

例如：blue，yellow，grey，lavender

根据喜好搭配

pch=节点的形状（从0~25，每个数字对应不同形状）

type参数用来控制所生成散点图的类型，有如下几个选项：

type=“p”表示绘制单独的点

type=“l”表示绘制点连成的折线

type=“b”表示有线连接的点

type=“o”表示将点绘在线上

type=“h”表示从点到x轴的垂线

type=“s”表示阶梯式图

type=“n”表示不绘制图

plot(x,main = "1900-1998年全球7级以上地震发生次数",xlab = "时间",ylab="7级以上地震发生次数",xlim = c(1900,1998),ylim=c(0,50),type="o",pch=24,col=5
)

04 读取包

包需要提前下载才可以读取

library(fBasics)     #加载包
library(fUnitRoots)
library(forecast)

下载包的流程：

install.packages(fBasics)

运行这个代码就可以了，不同包更改括号里的内容即可

05 单位根检验，又称平稳性检验

for(i in 0:1) print(adfTest(x,#adf检验，单位根检验lag=i,type="nc",title = "7级以上地震发生次数单位根检验"))
for(i in 0:3) print(adfTest(x,lag=i,type="c",title = "7级以上地震发生次数单位根检验"))
for(i in 0:3) print(adfTest(x,lag=i,type="ct",title = "7级以上地震发生次数单位根检验"))
adfTest(x)

代码运行结果：

> for(i in 0:1) print(
+ adfTest(x,#adf检验，单位根检验
+ lag=i,
+ type="nc",
+ title = "7级以上地震发生次数单位根检验"))

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 0
STATISTIC:
Dickey-Fuller: -1.5845
P VALUE:
0.1081

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -1.0451
P VALUE:
0.28

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

> for(i in 0:3) print(
+ adfTest(x,
+ lag=i,
+ type="c",
+ title = "7级以上地震发生次数单位根检验"))

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 0
STATISTIC:
Dickey-Fuller: -5.354
P VALUE:
0.01

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -3.9179
P VALUE:
0.01

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 2
STATISTIC:
Dickey-Fuller: -3.1832
P VALUE:
0.02451

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 3
STATISTIC:
Dickey-Fuller: -2.9558
P VALUE:
0.04427

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Warning messages:
1: In adfTest(x, lag = i, type = "c", title = "7级以上地震发生次数单位根检验") :
p-value smaller than printed p-value
2: In adfTest(x, lag = i, type = "c", title = "7级以上地震发生次数单位根检验") :
p-value smaller than printed p-value
> for(i in 0:3) print(
+ adfTest(x,
+ lag=i,
+ type="ct",
+ title = "7级以上地震发生次数单位根检验"))

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 0
STATISTIC:
Dickey-Fuller: -5.5514
P VALUE:
0.01

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -4.1418
P VALUE:
0.01

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 2
STATISTIC:
Dickey-Fuller: -3.5099
P VALUE:
0.04484

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Title:
7级以上地震发生次数单位根检验

Test Results:
PARAMETER:
Lag Order: 3
STATISTIC:
Dickey-Fuller: -3.3716
P VALUE:
0.06336

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

Warning messages:
1: In adfTest(x, lag = i, type = "ct", title = "7级以上地震发生次数单位根检验") :
p-value smaller than printed p-value
2: In adfTest(x, lag = i, type = "ct", title = "7级以上地震发生次数单位根检验") :
p-value smaller than printed p-value
> adfTest(x)

Title:
Augmented Dickey-Fuller Test

Test Results:
PARAMETER:
Lag Order: 1
STATISTIC:
Dickey-Fuller: -1.0451
P VALUE:
0.28

Description:
Mon Jan 9 15:40:36 2023 by user: Lenovo

从代码运行结果里找PVALUE的数值，数值大于0.05就是拒绝原假设，单位根不存在，序列平稳。

06 自相关图和偏相关图

par(mfrow=c(2,1))#分画布2行1列
acf(x)    #自相关
pacf(x)   #偏自相关

acf自相关（上） pacf偏自相关（下）

观察两图，发现自相关拖尾，偏自相关1阶截尾，预测模型为AR模型

07 白噪声检验或纯随机性检验

for(k in 1:2) print(      #白噪声检验,分析前后各一次Box.test(x,lag=6*k,            #type和l等价type="Ljung-Box"))   #服从卡方分布
#x-squared为卡方分布值，查表的值
#df为卡方检验自由度
#p-value 为p值
#若想显著性检验，检验目标为残差，验证残差是否为白噪声，为白噪声好

代码运行结果：

观察p-value的值，两值均小于0.05，拒绝原假设，序列非白噪声，序列有继续研究的价值。

08 模型拟合

x.fit<-arima(x,order=c(1,0,0),m="ML");x.fit#平稳序列模型

运行结果：

图1

根据自相关拖尾，偏相关1阶截尾，模拟AR(1)

x.fit适合x的一个模型

ML为最大似然估计法，CSS条件最小估计法，ML-CSS两种混合法

中间位置0固定 (差分)，若是MA(1)则c(0,0,1)

arima(x=序列x,)

Coefficients系数

ar1，模型AR系数

intercept截距miu

Xt的一个公式

s.e.估计量的标准差

sigma^2 estimated as 某某的估计是

log likelihood 极大似然对数值

aic最小信息准则，越小越好（多个模型间的比较）

09 显著性检验

re<-x.fit$residuals     #残差re
for(k in 1:2) print(      Box.test(re,       #残差lag=6*k,            type="Ljung-Box"))
qt(0.975,97)
t1=0.5432/0.0840;t1
t2=19.8911/1.3179;t2  #负数看左尾，正数看右尾

t1,t2的值可以对应上图1，不同模型对应输入，是一个2倍法则

代码运行结果：

观察p-value的值，两值均大于0.05，接受原假设，残差序列为白噪声序列，且t1,t2均大于1.984723,（2倍法则）说明拟合模型有效。

10 预测未来趋势

forecast为预测函数

h=_为预测的年数，h=5就是预测未来五年

xfore<-forecast(x.fit,h=5);xfore
#forecast预测
plot(xfore,col=2)
#出现的图形，蓝色的线是预测值，线两侧阴影为85%，95%的置信区间

代码运行结果：

未来五年置信区间

阴影部分为预测区间

结语：

以上关于时间序列图的绘制即检验、预测就完成了。

感谢大家阅读，如果方便的话给个赞同，谢谢！！！