ARIMA模型也可以对具有季节效应的序列建模。根据季节效应提取的方式不同，又分为ARIMA加法模型和ARIMA乘法模型。
        ABIMA加法模型是指序列中季节效应和其他效应之间是加法关系，即
                                                     
         这时，各种效应信息的提取都非常容易。通常简单的周期步长差分（即以周期为步长的步差分）即可将序列中的季节信息提取充分，简单的低阶差分（一般不超过三阶）即可将趋势信息提取充分，提取完季节信息和趋势信息之后的残差序列就是一个平稳序列，可以用ARMA模型拟合。          所以季节加法模型实际上就是通过趋势差分、季节差分将序列转化为平稳序列，再对其进行拟合。它的模型结构通常如下：
                                                   
式中，
        S为周期步长，d为提取趋势信息所用的差分阶数。
       为白噪声序列，且，。
        ，为q阶移动平均系数多项式。
        ，为p阶自回归系数多项式。
        该加法模型简记为，或。

以1962年—1991年德国工人季度失业率为例，用Rstudio软件进行模型拟合分析。数据如下，用时需保存为csv文件。

time x
1962/1/1 1.1
1962/4/1 0.5
1962/7/1 0.4
###### 0.7
1963/1/1 1.6
1963/4/1 0.6
1963/7/1 0.5
###### 0.7
1964/1/1 1.3
1964/4/1 0.6
1964/7/1 0.5
###### 0.7
1965/1/1 1.2
1965/4/1 0.5
1965/7/1 0.4
###### 0.6
1966/1/1 0.9
1966/4/1 0.5
1966/7/1 0.5
###### 1.1
1967/1/1 2.9
1967/4/1 2.1
1967/7/1 1.7
###### 2
1968/1/1 2.7
1968/4/1 1.3
1968/7/1 0.9
###### 1
1969/1/1 1.6
1969/4/1 0.6
1969/7/1 0.5
###### 0.7
1970/1/1 1.1
1970/4/1 0.5
1970/7/1 0.5
###### 0.6
1971/1/1 1.2
1971/4/1 0.7
1971/7/1 0.7
###### 1
1972/1/1 1.5
1972/4/1 1
1972/7/1 0.9
###### 1.1
1973/1/1 1.5
1973/4/1 1
1973/7/1 1
###### 1.6
1974/1/1 2.6
1974/4/1 2.1
1974/7/1 2.3
###### 3.6
1975/1/1 5
1975/4/1 4.5
1975/7/1 4.5
###### 4.9
1976/1/1 5.7
1976/4/1 4.3
1976/7/1 4
###### 4.4
1977/1/1 5.2
1977/4/1 4.3
1977/7/1 4.2
###### 4.5
1978/1/1 5.2
1978/4/1 4.1
1978/7/1 3.9
###### 4.1
1979/1/1 4.8
1979/4/1 3.5
1979/7/1 3.4
###### 3.5
1980/1/1 4.2
1980/4/1 3.4
1980/7/1 3.6
###### 4.3
1981/1/1 5.5
1981/4/1 4.8
1981/7/1 5.4
###### 6.5
1982/1/1 8
1982/4/1 7
1982/7/1 7.4
###### 8.5
1983/1/1 10.1
1983/4/1 8.9
1983/7/1 8.8
###### 9
1984/1/1 10
1984/4/1 8.7
1984/7/1 8.8
###### 8.9
1985/1/1 10.4
1985/4/1 8.9
1985/7/1 8.9
###### 9
1986/1/1 10.2
1986/4/1 8.6
1986/7/1 8.4
###### 8.4
1987/1/1 9.9
1987/4/1 8.5
1987/7/1 8.6
###### 8.7
1988/1/1 9.8
1988/4/1 8.6
1988/7/1 8.4
###### 8.2
1989/1/1 8.8
1989/4/1 7.6
1989/7/1 7.5
###### 7.6
1990/1/1 8.1
1990/4/1 7.1
1990/7/1 6.9
###### 6.6
1991/1/1 6.8
1991/4/1 6
1991/7/1 6.2
###### 6.2

拟合过程如下：

1、绘制观察值序列时序图

由于观察值是以季度为单位进行记录，也就是一年中记录4次数据，所以生成时间序列时frequency为4。

m<-read.csv("C:\\Users\\86191\\Desktop\\A1_20.csv",sep=',',header = T)
x<-ts(m$x,start = c(1962,1),frequency = 4)#绘制时序图
plot(x,main = "1962年-1991年德国工人季度失业率",xlab = "年份",ylab="gov_cons",xlim = c(1962,1992),ylim=c(0,11),type="o",pch=20,
)

绘制时序图如下：

可以看到，该序列既含有长期趋势又含有以年为周期的季节效应。不是平稳序列，接下来需要对其提取趋势，季节效应。

        2.差分平稳化

将第1步中的序列记作原序列，则对原序列作1阶差分消除(提取)趋势，再作4步差分消除季节效应的影响。

#绘制一阶差分再以周期时长数目为k的k步差分的时序图
x_diff_4<-diff(diff(x),4)#一阶差分再4步差分
plot(x_diff_4,main = "1962年-1991年德国工人季度失业率一阶差分再4步差分",type="o",pch=20,
)

差分后的序列的时序图如下：

时序图显示，差分后的序列没有明显的趋势或周期，随机波动比较平稳。但这仅是图示法，具有一定的主观性，需要进行ADF检验，白噪声检验判断序列是否具有分析价值。

从时序图可以看到，序列大致在0附近来回波动，所以选择nc类型(无漂移项自回归结构)的ADF检验。

#对一阶差分再以周期时长数目为k的k步差分后的序列进行平稳性检验
library(fBasics)
library(fUnitRoots)
for(i in 0:3) print(adfTest(x_diff_4,lag=i,type="nc",title = "1962年-1991年德国工人季度失业率单位根检验"))

结果如下：

Title:
 1962年-1991年德国工人季度失业率单位根检验

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

Description:
 Sun Nov 13 00:24:28 2022 by user: 86191

Title:
 1962年-1991年德国工人季度失业率单位根检验

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

Description:
 Sun Nov 13 00:24:28 2022 by user: 86191

Title:
 1962年-1991年德国工人季度失业率单位根检验

Test Results:
  PARAMETER:
    Lag Order: 2
  STATISTIC:
    Dickey-Fuller: -4.8907
  P VALUE:
    0.01

Description:
 Sun Nov 13 00:24:28 2022 by user: 86191

Title:
 1962年-1991年德国工人季度失业率单位根检验

Test Results:
  PARAMETER:
    Lag Order: 3
  STATISTIC:
    Dickey-Fuller: -6.2638
  P VALUE:
    0.01

Description:
 Sun Nov 13 00:24:28 2022 by user: 86191

Warning messages:
1: In adfTest(x_diff_4, lag = i, type = "nc", title = "1962年-1991年德国工人季度失业率单位根检验") :
  p-value smaller than printed p-value
2: In adfTest(x_diff_4, lag = i, type = "nc", title = "1962年-1991年德国工人季度失业率单位根检验") :
  p-value smaller than printed p-value
3: In adfTest(x_diff_4, lag = i, type = "nc", title = "1962年-1991年德国工人季度失业率单位根检验") :
  p-value smaller than printed p-value
4: In adfTest(x_diff_4, lag = i, type = "nc", title = "1962年-1991年德国工人季度失业率单位根检验") :
  p-value smaller than printed p-value

结果显示，0-3阶延迟的ADF检验，其p值都显示为0.01，但提示信息告诉我们，真实的p值小于0.01，因此在显著性水平为0.01下，没有充分理由接受ADF检验的原假设：序列非平稳；从而认为差分后的序列是平稳的。接下来进行白噪声检验。

#对一阶差分再以周期时长数目为k的k步差分后的序列进行白噪声检验
for(k in 1:2) print(Box.test(x_diff_4,lag=6*k,type="Ljung-Box"))

结果如下：

Box-Ljung test

data:  x_diff_4
X-squared = 43.837, df = 6, p-value = 7.964e-08

Box-Ljung test

data:  x_diff_4
X-squared = 51.708, df = 12, p-value = 6.982e-07

检验结果显示，延迟6阶与延迟12阶的白噪声检验的p值都远小于0.001，没有充分理由接受白噪声检验的原假设：序列是白噪声序列；从而认为差分后的序列是非白噪声的。                                  也就是说，差分后的序列是平稳非白噪声序列，具有分析价值，接下来对差分后的序列进一步拟合ARMA模型。

      3、模型拟合

        考察差分后序列的自相关图和偏自相关图的性质，给拟合模型定阶。

acf(x_diff_4)#绘制一阶差分再以周期时长数目为k的k步差分自相关图
pacf(x_diff_4)#绘制一阶差分再以周期时长数目为k的k步差分偏相关图

自相关图如下：        偏自相关图如下：
        自相关图显示出明显的下滑轨迹，这是典型的拖尾属性。偏自相关图除了1阶和4阶偏自相关系数显著大于2倍标准差，其他阶数的偏自相关系数基本都在2倍标准差范围内波动。所以尝试拟合疏系数模型AR(1,4)。由于第1步中对原序列进行了差分，所以实际上是拟合疏系数的季节加法模型ARIMA((1,4),(1,4),0)，或等价表达为ARIMA(1,1,0)(0,1,0)。使用条件最小二乘法估计模型参数。

#模型拟合-加法季节模型
x_fit<-arima(x,order = c(4,1,0),seasonal=list(order=c(0,1,0),period=4),method = "CSS",transform.pars = F,fixed = c(NA,0,0,NA))#实际拟合模型
#include.mean = T：拟合常数项
#include.mean = F：不拟合常数项
#transform.pars：参数估计是否由系统完成
#transform.pars = T：系统根据order设置的模型阶数自动完成参数估计
#transform.pars = F：需要拟合疏系数模型，需要认为干预
#fixed：对疏系数模型指定疏系数的位置
#seasonal:指定季节模型的阶数与季节周期，命令格式为seasonal=list(c(p,d,q),period=周期),d代表以周期为步长的步差分次数
#乘法模型中p,q不全为零,加法模型中p=0,q=0
x_fit

结果如下：

Call:
arima(x = x, order = c(4, 1, 0), seasonal = list(order = c(0, 1, 0), period = 4), 
    transform.pars = F, fixed = c(NA, 0, 0, NA), method = "CSS")

Coefficients:
             ar1  ar2  ar3      ar4
        0.4507    0    0  -0.2816
s.e.  0.0804    0    0   0.0808

sigma^2 estimated as 0.09377:  part log likelihood = -27.08

得到该模型的口径如下：                  

4、模型检验

对差分后的序列的残差进行白噪声检验，以判断拟合的模型是否有效，即序列信息是否提取充分。

#模型的显著性检验(有效性)-残差序列随机性检验
for(k in 1:2) print(Box.test(x_fit$residuals,lag=6*k,type="Ljung-Box"))

结果如下：

Box-Ljung test

data:  x_fit$residuals
X-squared = 1.7605, df = 6, p-value = 0.9404

Box-Ljung test

data:  x_fit$residuals
X-squared = 9.1286, df = 12, p-value = 0.6919

检验结果显示，延迟6阶与延迟12阶的白噪声检验的p值都远大于0.05的显著性水平，从而认为残差序列为白噪声序列，信息提取充分。接下来对模型的参数进行显著性检验。

#参数显著性检验，即检验参数是否显著非零,拒绝原假设则显著非零
qt(0.975,118)#检验统计量临界值
#ar1系数检验
t1<-0.4507/0.0804
t1
pt(t1,df=118,lower.tail = F)#ar4检验
t2<-0.2816/0.0808
t2
pt(t2,df=188,lower.tail = T)

检验结果如下：

> #参数显著性检验，即检验参数是否显著非零,拒绝原假设则显著非零
> qt(0.975,118)#检验统计量临界值
[1] 1.980272
> #ar1系数检验
> t1<-0.4507/0.0804
> t1
[1] 5.605721
> pt(t1,df=118,lower.tail = F)
[1] 6.921003e-08
> 
> #ar4检验
> t2<-0.2816/0.0808
> t2
[1] 3.485149
> pt(t2,df=188,lower.tail = T)
[1] 0.999694

结果显示，两参数的检验统计量值的绝对值均大于临界值(n为序列长度，m为待估参数个数)，即模型参数都显著非零，而第二个参数的检验p值为0.999694是因为R语言中的定义，而并非其p值真是0.999694，可根据其检验统计量的值与临界值判断。
        综述，结果显示该模型拟合良好，对序列相关信息的提取充分。
        将序列拟合值和序列观察值联合作图，并基于该模型预测1992-2001年德国工人失业率。

#预测
library(forecast)
x.fore<-forecast(x_fit,h=20,level = c(95))#h为预测期数,level为置信度,预测五年
x.fore
plot(x.fore)#详细预测图-星号是观察值序列，蓝色是两倍95%置信线，实线是拟合值序列
L1<-x.fore$fitted-1.96*sqrt(x_fit$sigma2)
U1<-x.fore$fitted+1.96*sqrt(x_fit$sigma2)
L2<-ts(x.fore$lower[,1],start=1992)#start=预测期数起始年限
U2<-ts(x.fore$upper[,1],start=1992)#start=预测期数起始年限
c1<-min(x,L1,L2)
c2<-max(x,L2,U2)
plot(x,type="p",pch=8,xlim=c(1962,1995),ylim=c(c1,c2))#xlim=c(数据序列起始时间,预测期数终止时间)
lines(x.fore$fitted,col=2,lwd=2)#此时拟合值与观测值同图
lines(x.fore$mean,col=2,lwd=2)
lines(L1,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(L1,col=4,lty=2)
lines(L2,col=4,lty=2)
lines(U2,col=4,lty=2)

1992-2001年德国工人的失业率预测结果如下：

Point Forecast             Lo 95        Hi 95
1992 Q1       6.619702  6.01950988  7.219894
1992 Q2       5.862413  4.80487129  6.919955
1992 Q3       5.969028  4.51858464  7.419471
1992 Q4       5.842457  4.05323318  7.631681
1993 Q1       6.143241  3.80612560  8.480356
1993 Q2       5.320323  2.46925775  8.171388
1993 Q3       5.423652  2.12609740  8.721208
1993 Q4       5.331242  1.64292433  9.019560
1994 Q1       5.680910  1.42577545  9.936045
1994 Q2       4.898508  0.06288166  9.734133
1994 Q3       5.021024 -0.36391902 10.405967
1994 Q4       4.927643 -0.96927958 10.824565
1995 Q1       5.263108 -1.30239979 11.828616
1995 Q2       4.462895 -2.77857660 11.704366
1995 Q3       4.571980 -3.30930009 12.453260
1995 Q4       4.472818 -4.00579982 12.951436
1996 Q1       4.809677 -4.40046864 14.019823
1996 Q2       4.015107 -5.93377754 13.963992
1996 Q3       4.130519 -6.52681951 14.787857
1996 Q4       4.035836 -7.29318279 15.364855

预测图如下：        通过图示也可以直观地看出该模型对序列的拟合效果良好。

应用时间序列分析——有季节效应的非平稳序列分析-ARIMA加法模型-R语言相关推荐

  1. 应用时间序列分析——有季节效应的非平稳序列分析-R语言

    因素分解理论: 1919年统计学家沃伦·珀森斯(Warren Persons)在他的论文<商业环境的指标>中首次提出了确定性因素分解(time series decomposition)思 ...

  2. 应用时间序列分析--有季节效应的非平稳序列分析-ARIMA乘法模型-R语言

    在上一篇文章(http://t.csdn.cn/Lp2Nj)中,我们用到的数据是一个既含有季节效应又含有长期趋势效应的简单序列,为什么说它简单呢,是因为这种序列的季节效应.趋势效应和随机波动彼此之间很 ...

  3. 有季节效应的非平稳序列分析之因素分解

    目录 6.1 因素分解理论 6.2因素分解模型 6.2.1因素分解模型的选择 6.2.2趋势效应的提取 简单中心移动平均的良好属性 R语言中,使用filter函数可以做简单移动平均 6.2.3 季节效 ...

  4. 无季节效应的非平稳序列分析（一)

    目录 Cramer分解定理(1961年提出) 差分 R语言函数 diff 例题: 过差分: 小结 Cramer分解定理(1961年提出) 任何一个时间序列 都可以分解为两部分的叠加:其中一部分是由多项 ...

  5. 时间序列、时间序列分析、时间序列效应分解、平稳时间序列、AIC和BIC

    时间序列.时间序列分析.时间序列效应分解.平稳时间序列.AIC和BIC 目录 时间序列.时间序列分析.时间序列效应分解.平稳时间序列.AIC和BIC

  6. 使用连续小波变换(cwt)对时间轴上频率发生变化的非平稳信号进行处理

    使用连续小波变换对时间轴上频率发生变化的非平稳信号进行处理 傅里叶变换对时间轴上的平稳信号处理有着很大的优势,但是对非平稳信号就无能为力了. 小波变换是基于傅里叶变换产生的可以用于处理非平稳信号的信号 ...

  7. R语言逻辑运算符（Logical Operators，大于、小于、等于、不等于、与或非、是否为真）、R语言逻辑运算符（Logical Operators）实战示例

    R语言逻辑运算符(Logical Operators,大于.小于.等于.不等于.与或非.是否为真.>.<.!=.==.&.|.!&&.||).R语言逻辑运算符(Lo ...

  8. R语言ARIMA集成模型预测时间序列分析

    全文链接:http://tecdat.cn/?p=18493 本文我们使用4个时间序列模型对每周的温度序列建模.第一个是通过auto.arima获得的,然后两个是SARIMA模型,最后一个是Buys- ...

  9. 决策树模型回归可视化分析_【时间序列分析】在论文中用向量自回归（VAR）模型时应注意哪些问题？...

    在论文的写作中,向量自回归(VAR)模型是经常用的一个模型,同时它也是多维时间序列模型的最核心内容之一. 首先要清楚,VAR模型主要是考察多个变量之间的动态互动关系,从而解释各种经济冲击对经济变量形成 ...

  10. R语言线性回归和时间序列分析北京房价影响因素可视化案例

    原文链接:http://tecdat.cn/?p=21467 目的 房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子 人口众多 独生子女政策:如何影响房子的几何 ...

最新文章

  1. pcb 假八层_八层板的成本却只六层的性能，PCB设计遭遇假八层该怎么办？
  2. java-数组排序--冒泡排序、鸡尾酒排序、地精排序
  3. 大志非才不就，大才非学不成—博文资源汇总
  4. oracle获取当前日期的前一天，判断两个日期是否相等
  5. 值重新赋值_JavaScript-赋值运算符
  6. JavaScript匿名函数以及在循环中的匿名函数
  7. mysql体系结构和存储引擎
  8. java读取同包文件_Java实现从jar包中读取指定文件的方法
  9. 服务器怎么用u盘传文件进去_云服务器无法使用USB设备？一文读懂如何在云服务器上使用优盘...
  10. route-policy和ACL组合时permit和deny的作用
  11. 字符串中最长不重合子串长度
  12. mysql管理命令 java_Mysql常用命令及操作
  13. html 超链接嵌套,嵌套的超链接区域,HTML源中没有嵌套的链接元素
  14. Spring的两种定时器
  15. Packet Tracer 思科模拟器入门教程 之十 路由器单臂路由配置
  16. 介绍一种养眼的IDEA/Pycharm/Clion配色方案
  17. SM2证书的鉴定方法——续上文国密算法
  18. 干货分享：数据分析面试一般都问啥？
  19. properties文件的用法；utils类封装数据库连接、资源关闭、通用的增删改，以及两种通用的查询方法
  20. 大带宽服务器对于网站速度的影响有多大？

热门文章

  1. NIOS 2 软核中EPCS配置芯片的存储操作
  2. 通过python实现同步修改本地电脑时间
  3. matlab中在xls单元格中填充颜色,!Excel中如何根据某一列的值对整行进行颜色填充？...
  4. 探索图片透明度叠加方式
  5. Ddos攻击怎么防护？DDOS八大防御策略
  6. 蓝光三维扫描仪用于无人机结构外壳虚拟装配
  7. 软件测试培训到底值不值得参加？
  8. 318公路是中国最长最美国道，沿途风景绝佳，进藏必去的最佳路线
  9. word 大纲视图的生成
  10. 鸿蒙大陆v2.8正式版,鸿蒙大陆正式版地图下载-鸿蒙大陆下载 V2.22--pc6下载站