精华内容
下载资源
问答
  • R语言时间序列分析常用步骤
    千次阅读
    2021-05-12 17:40:52

     常用步骤代码:

    #载入数据
    d <- WWWusage
    x <- ts(d,start = 1) # 其实WWWusage本身就是时间序列格式的,这里是为了提醒大家记得ts
    plot(x)
    
    #差分并观察
    x.dif <- diff(x,1,2)
    plot(x.dif)
    #差分序列ADF检验
    library(tseries)
    adf.test(x.dif) # p越小越好
    #确定差分阶数d
    
    #序列纯随机性检验
    for(i in 1:2)print(Box.test(x.dif, lag=6*i))
    #p大,拒绝序列为白噪声的假设,说明序列有规律可提取,有研究价值,可继续建模
    
    #判断p、q
    acf(x.dif)
    pacf(x.dif)
    #根据auto.arima辅助判断
    library(forecast)
    auto.arima(x.dif) # AIC 实际值越小越好
    x.fit <- arima(x.dif,order = c(2,0,0));x.fit
    #多对比几次结果
    
    #AIC是赤池统计量,当增加了一个变量后,如果这个变量对模型有足够的解释力,就会使赤池统计量变小
    #AIC的值为负的时候也是越小越好,也就说是实际值而不是绝对值越小越好,比如-2与-1.9,就选-2那个
    #AIC和BIC是似然函数和乘法函数的线性组合,不同之处在于惩罚力度不同。都是越小越好
    #loglik是对数似然,模型比较时越大越好
    
    #残差白噪声检测,考察拟合效果,p越大,说明残差纯随机性越显著,信息充分提取
    for(i in 1:2) print(Box.test(x.fit$residuals,lag = 6*i))
    
    #通过检测的模型进行预测
    library(zoo)
    library(forecast)
    x.fore <- forecast(x.fit,h = 10) #预测未来10期
    plot(x.fore)
    x.fore

    模型阶数判定:

    以datasets包自带WWWusage数据为例,2阶差分后拟合模型,AIC实际值(不是绝对值)越小越好,实验如下:

    win.graph();par(mfrow=c(1,2));acf(x.dif);pacf(x.dif)

    结合acf图和pacf图,结果似乎是p=2,q=2最佳,但对比AIC值却不是如此:

     结果表明对x <- ts(WWWusage)来说,arima(2,2,0)的拟合效果最好,AIC值为511.46。

    为什么偏自相关图判断的结果不是拟合效果最好的呢?希望哪位老师同学可以指点我,谢谢!

    我的微信:1961812312


    更新:

    之前没有对p,q作为组合代入模型寻找最小aic值,

    下一篇文章中,分析了36种组合的结果,找到aic值最小的结果为p=5,q=5,aic值为509.8135

    但还是与自相关图和偏自相关图结果给出的建议不一样,这一点还是希望有老师同学指点一下!


    参考文献:

    https://people.duke.edu/~rnau/411arim.htm

    https://bbs.pinggu.org/forum.php?mod=viewthread&action=printable&tid=226732

    https://bbs.pinggu.org/thread-3226868-1-1.html

    https://www.zybuluo.com/evilking/note/851533

    wikimirror

    感谢上帝的恩典。

    另:根据ARIMA模型得到公式结果,可以参考wiki和Google的镜像资料代入得出。


    更新:

    R语言时间序列分析-根据aic值选择arima模型


    更新2:

    对于之前的问题,当(偏)自相关图给出的模型不是AIC值最优的模型时,应该选择哪一个?答案如下:

    从ACF和PACF图上定性判断有时不准确,建议以量化的标准(AIC值)为准。不放心也可以再多比较一下其他指标,包括:BIC值,两个模型的残差随机性检验,最佳子集等

    ——网友何劼

    感谢网友@没有那把剑、@ 何劼

    更多相关内容
  • R语言时间序列分析

    千次阅读 多人点赞 2020-05-23 00:29:48
    当拿到一个时间序列的时候,首先分析时间序列的类型,不同类型的序列有不同的处理方式。本文包含以下几个部分: 1、时间序列数据准备 2、时间序列平稳性检验 3、拟合ARIMA模型 4、ARIMA模型的检验诊断 ...

    R语言时间序列分析


    时间序列白噪声序列终止分析,无信息可提取
    非白噪声序列平稳序列AR/MA/ARMA
    非平稳序列ARIMA,实际应用中最常见

    当拿到一个时间序列的时候,首先分析该时间序列的类型情况,不同类型的序列有不同的处理方式。本文包含以下几个部分:

    1、时间序列数据准备
    2、时间序列平稳性检验
    3、拟合ARIMA模型
    4、ARIMA模型的检验诊断
    5、用ARIMA模型进行预测

    一、时间序列数据准备

    > stodata <- read.csv("Stk_Day.csv")
    > head(stodata)
          代码       时间 开盘价 最高价 最低价 收盘价 成交量.股. 成交额.元.
    1 SH600000 1999-11-10  29.50  29.80  27.00  27.75  174085000 4859102208
    2 SH600000 1999-11-11  27.58  28.38  27.53  27.71   29403400  821582016
    3 SH600000 1999-11-12  27.86  28.30  27.77  28.05   15007900  421591008
    4 SH600000 1999-11-15  28.20  28.25  27.70  27.75   11921000  332952000
    5 SH600000 1999-11-16  27.88  27.97  26.48  26.55   23223100  628908032
    6 SH600000 1999-11-17  26.50  27.18  26.37  27.18   10052500  268995008
    > library(xts)
    > x1 <- stodata[,3]
    > x2 <- stodata[,7]
    > dt <- as.Date(stodata[,2]) #dt必须是时间格式
    > tsdata1 <- xts(x1,dt)
    > tsdata2 <- xts(x2,dt) #建立时间序列类型
    > par(mfrow=c(2,1)) #画图时两行一列的排列
    > plot(tsdata1,col="black")
    > plot(tsdata2,col="red")

     

    二、时间序列平稳性检验

    (1)时序图主观检验

    平稳时间序列的均值和方差都为常数。根据这一性质,平稳序列的时序图显示该序列值始终在一个常数附近随机波动,而且波动的范围有界;如果有明显的趋势性或周期性那它通常不是平稳序列。观察tsdata1的时间序列图,明显有趋势性,所以这个时间序列不是平稳序列(一般非平稳序列经过一定的差分计算会呈现出一定的平稳性)。

    (2)绘制自相关函数图检验

    > acf(tsdata1)

    平稳序列具有短期相关性,这个性质表明对平稳序列而言通常只有近期的序列值对现时值的影响比较明显,间隔越远的过去值对现时值的影响越小。随着延迟期数k的增加,平稳序列的自相关系数ρ(k)会比较快的衰减趋向于零,并在零附近随机波动(平稳时序的自相关系数具有拖尾性,即模型的自相关系数ρ(k)呈指数的速度衰减,始终有非零取值,不会在k大于某个常数之后就恒等于零)。而非平稳序列的自相关系数衰减的速度比较慢,这就是利用自相关图进行平稳性检验的标准。

    可以通过自相关和偏自相关图来判断时间序列是否平稳。平稳的序列的自相关图和偏相关图不是拖尾就是截尾。截尾是指在某阶之后,系数都为 0;拖尾是指有一个衰减的趋势,但是不都为0。

    ACF图形状序列模型
    不衰减到0不平稳
    固定点上有大于0的值(周期性)周期性
    所有点都是0完全随机(白噪声)
    指数衰减到0(拖尾)AR(p)
    正负交替衰减到0(拖尾)AR(p)
    前几个大于0,后面都为0(截尾)MA(q)

    图中的水平虚线是相关系数是否显著的分界线,在虚线之上的是显著的。显著的自相关是考虑应用自回归移动平均模型ARIMA来对时间序列建模的一个证据。可以看出,tsdata1不是平稳序列,但是是自相关的。

    > d_tsdata1 <- diff(tsdata1) #diff函数计算差分,输出还是一个时间序列格式
    > acf(d_tsdata1[-1]) #-1去掉差分后的第一项NA

    tsdata1差分后的时间序列d_tsdata1明显具有“拖尾性”,可见tsdata1的差分序列是平稳的。这也是ARIMA模型原理的关键所在之处。

    (3)检验时间序列的自相关

    以上两种“看图”的方式都是有人为主观参与的检验方式,其实可以用严格的显著性检验来得到结论的。Box.test(ts)函数输出一个p值,一般情况下,如果p值小于0.05,说明时间序列有显著的自相关;而p值大于0.05,则没有证据证明自相关存在。本例中,p值<<0.05证明时间序列tsdata1是自相关的。

    > Box.test(tsdata1)
    
    	Box-Pierce test
    
    data:  tsdata1
    X-squared = 3670.1, df = 1, p-value < 2.2e-16

    最后总结:自相关的序列一定是 [平稳序列] 或者 [非平稳序列],肯定不是 [白噪声序列] ,通过自相关检验的序列一定可以用ARIMA模型拟合。

    三、拟合ARIMA模型

    > library(forecast)
    Registered S3 method overwritten by 'quantmod':
      method            from
      as.zoo.data.frame zoo 
    > auto.arima(tsdata1) #ARIMA模型自动拟合
    Series: tsdata1 
    ARIMA(5,1,1) 
    
    Coefficients:
             ar1     ar2     ar3     ar4      ar5      ma1
          0.5378  0.0076  0.0145  0.0415  -0.0980  -0.5958
    s.e.  0.0830  0.0192  0.0188  0.0187   0.0165   0.0823
    
    sigma^2 estimated as 0.3692:  log likelihood=-3390.94
    AIC=6795.88   AICc=6795.91   BIC=6839.37
    > m <- arima(tsdata1,order=c(5,1,1)) #如果提前知道最优阶数,也可以用直接用ARIMA函数,两种方式等价
    > confint(m) #展示ARIMA模型系数的置信区间
               2.5 %      97.5 %
    ar1  0.375048689  0.70050206
    ar2 -0.030063447  0.04527976
    ar3 -0.022297674  0.05135013
    ar4  0.004946544  0.07812118
    ar5 -0.130321386 -0.06562303
    ma1 -0.757130275 -0.43451464
    > m <- arima(tsdata1,order=c(5,1,1),fixed=c(NA,0,0,NA,NA,NA)) #用fixed参数提出模型中不显著的系数
    Warning message:
    In arima(tsdata1, order = c(5, 1, 1), fixed = c(NA, 0, 0, NA, NA,  :
      一些AR参数是固定的:把transform.pars设成FALSE

    模型的阶数由3个整数来表示,(p,d,q),这里的p是自回归系数的个数,d是差分的阶数,q是移动平均系数的个数。当建立arima模型时,一般不知道什么是合适的阶数。一般会用函数auto.arima(),让它来为我们选择合适的模型阶数,而不是手动繁琐地寻找p,d,q的最优组合。

    ARIMA模型有一个令人头痛的问题:并不是所有的系数都是显著,也就是回归的系数可能等于0。当存在不显著的系数的时候就需要进行剔除,函数arima()含有一个向量参数fixed,可以对ARIMA模型中不显著的系数进行剔除。在本例中,ar2,ar3的置信区间包含0值,因此可以剔除。这里函数输出的警告信息不用管。

    四、ARIMA模型的检验诊断

    > tsdiag(m)

    函数tsdiag()给出一组(3个)图形。通过对这三个图形进行分析可对ARIMA模型进行诊断:

    (1)标准化残差图。一个好的ARIMA模型的标准化残差没有波动聚集性,本例中具有一定的聚集性,不过仍可用于预测分析。

    (2)残差自相关图。一个好的ARIMA模型的自相关函数没有显著的自相关性,本例中的ACF图表现很好,没有显著的自相关性。

    (3)Ljung-Box统计量p值图。如果Ljung-Box统计量的p值都比较大,表明残差没有任何模式,模型已经提取了数据中的所有信息,剩余的仅仅是“噪声”。本例是一个表现很好的例子。

    五、用ARIMA模型进行预测

    > p <- predict(m,n.ahead=50)
    > p
    $pred
    Time Series:
    Start = 3688 
    End = 3737 
    Frequency = 1 
     [1] 16.47158 16.63245 16.68664 16.60102 16.67366 16.71499 16.72320 16.71790
     [9] 16.72714 16.72682 16.72301 16.71997 16.71938 16.71816 16.71737 16.71718
    [17] 16.71736 16.71745 16.71757 16.71771 16.71780 16.71784 16.71785 16.71786
    [25] 16.71785 16.71784 16.71783 16.71782 16.71782 16.71782 16.71782 16.71782
    [33] 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782
    [41] 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782 16.71782
    [49] 16.71782 16.71782
    
    $se
    Time Series:
    Start = 3688 
    End = 3737 
    Frequency = 1 
     [1] 0.6072344 0.8367099 1.0064440 1.1473611 1.2840543 1.3885973 1.4779744
     [8] 1.5589136 1.6348650 1.7055611 1.7734991 1.8394046 1.9034873 1.9657107
    [15] 2.0262971 2.0852772 2.1426909 2.1986054 2.2531249 2.3063352 2.3583241
    [22] 2.4091757 2.4589676 2.5077667 2.5556329 2.6026193 2.6487736 2.6941385
    [29] 2.7387531 2.7826530 2.8258714 2.8684387 2.9103836 2.9517324 2.9925098
    [36] 3.0327389 3.0724413 3.1116372 3.1503454 3.1885837 3.2263689 3.2637167
    [43] 3.3006419 3.3371585 3.3732799 3.4090186 3.4443864 3.4793948 3.5140544
    [50] 3.5483755
    

    函数predict()会根据模型计算未来的观测值和他的标准误差。返回值是一个有两个元素的列表,一个元素pred是预测值,另一个元素se是标准误差。这里的pred和se都编码为时间序列格式,因此R会打印预测的数值,同时还有他们的时间序列属性,如果需要预测未来多个观测值,需要设置参数n.ahead,本例中向后预测50个观测值。需要注意的是,每向后预测一个,标准误差会变大,这是因为越往后预测就会变得越不确定。

    本例的预测效果不是太好,主要表达使用过程。


    tseries

    上面的时间序列分析的一些原理以及xts包,其实做时间序列分析还有许多包,下面记一下tseries包。

    一、数据准备与查看

    > library(tseries)
    
        ‘tseries’ version: 0.10-47
    
        ‘tseries’ is a package for time series analysis and
        computational finance.
    
        See ‘library(help="tseries")’ for details.
    > library(forecast)
    > air <- AirPassengers #关于ts数据格式的构建,可查询资料学习,本例直接应用一个给出的ts格式
    > air
         Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
    1949 112 118 132 129 121 135 148 148 136 119 104 118
    1950 115 126 141 135 125 149 170 170 158 133 114 140
    1951 145 150 178 163 172 178 199 199 184 162 146 166
    1952 171 180 193 181 183 218 230 242 209 191 172 194
    1953 196 196 236 235 229 243 264 272 237 211 180 201
    1954 204 188 235 227 234 264 302 293 259 229 203 229
    1955 242 233 267 269 270 315 364 347 312 274 237 278
    1956 284 277 317 313 318 374 413 405 355 306 271 306
    1957 315 301 356 348 355 422 465 467 404 347 305 336
    1958 340 318 362 348 363 435 491 505 404 359 310 337
    1959 360 342 406 396 420 472 548 559 463 407 362 405
    1960 417 391 419 461 472 535 622 606 508 461 390 432
    > plot(air)
    > tsdisplay(air) #直接使用tsdisplay来观察,它包含了时序图,以及acf、pacf两个相关图
    > air_train <- ts(as.vector(air[1:132]),frequency=12,start=c(1949,1)) #训练集
    > air_test <- ts(as.vector(air[133:144]),frequency=12,start=c(1960,1)) #测试集

    二、数据处理与检验

    > tsdisplay(air_train)
    #air_train存在一个向上的趋势,采用一阶差分消除
    > air_train1 <- diff(air_train,1)
    > tsdisplay(air_train1)
    #差分后的时间序列基本围绕在0波动,基本平稳
    > adf.test(air_train1) #单位根检验,p<0.05,序列air_train1为平稳序列
    
            Augmented Dickey-Fuller Test
    
    data:  air_train1
    Dickey-Fuller = -6.679, Lag order = 5, p-value = 0.01
    alternative hypothesis: stationary
    
    Warning message:
    In adf.test(air_train1) : p-value smaller than printed p-value 
    #差分后的序列acf图中显示存在一个12阶的季节性影响因素,然后通过做12阶差分图看一下是否可消除
    > tsdisplay(diff(air_train1,12))

    三、模型拟合与诊断

    > mo <- auto.arima(air_train)
    > mo
    Series: air_train 
    ARIMA(1,1,0)(0,1,0)[12] 
    
    Coefficients:
              ar1
          -0.2431
    s.e.   0.0894
    
    sigma^2 estimated as 109.8:  log likelihood=-447.95
    AIC=899.9   AICc=900.01   BIC=905.46
    > confint(mo) #查看模型系数,检验是否有需要剔除的系数
             2.5 %      97.5 %
    ar1 -0.4184353 -0.06782791
    > tsdiag(mo) #进行模型诊断,三个图表明是个很好的模型
    > qqnorm(mo$residuals)
    > qqline(mo$residuals) #用QQ图来检验残差的正态性,从而确定模型质量。如果所有点都落在45度线上表明模型表达非常好。
    > summary(mo)
    Series: air_train 
    ARIMA(1,1,0)(0,1,0)[12] 
    
    Coefficients:
              ar1
          -0.2431
    s.e.   0.0894
    
    sigma^2 estimated as 109.8:  log likelihood=-447.95
    AIC=899.9   AICc=900.01   BIC=905.46
    
    Training set error measures:
                       ME     RMSE      MAE       MPE     MAPE      MASE
    Training set 0.579486 9.907267 7.483159 0.1187348 2.880429 0.2457523
                       ACF1
    Training set 0.01227544
    

    四、模型预测与评估

    > pred <- forecast(mo,h=12,level=c(99.5))
    > pred
             Point Forecast  Lo 99.5  Hi 99.5
    Jan 1960       424.1099 394.6963 453.5234
    Feb 1960       407.0557 370.1672 443.9442
    Mar 1960       470.8257 426.8166 514.8349
    Apr 1960       460.8817 410.9544 510.8090
    May 1960       484.8681 429.6094 540.1268
    Jun 1960       536.8714 476.7621 596.9807
    Jul 1960       612.8706 548.2716 677.4695
    Aug 1960       623.8708 555.0751 692.6664
    Sep 1960       527.8707 455.1199 600.6215
    Oct 1960       471.8707 395.3690 548.3725
    Nov 1960       426.8707 346.7936 506.9479
    Dec 1960       469.8707 386.3711 553.3704
    > plot(pred)
    > lines(air_test,col="red") #将预测和实际值画在一个坐标图中,图中红线是实际值,蓝线是预测值


    推荐几篇很有用的文章,同时也是本文主要参考的文章:

    1.https://blog.csdn.net/gdyflxw/article/details/55509656

    2.https://zhuanlan.zhihu.com/p/27919755

    3.https://blog.csdn.net/LiuYuan_BJTU/article/details/67068404

    4.https://zhuanlan.zhihu.com/p/27453841

    展开全文
  • 时间序列分析代码.R

    2020-04-17 21:36:07
    时间序列分析R语言的相关代码,如AR模型,MA模型,ARMA模型,以及自相关检验,异方差检验,画自相关图,偏自相关图,函数定阶,模型系数的显著性检验,模型预测,输出预测图和拟合图等。
  • 利用R语言对化学浓度读数数据进行时间序列分析,建立了ARMA模型。附有全部代码以及相关数据集。
  • 自由选取一组数据(可以是R 自带的数据集、或者其它来源,鼓励选取一些有趣的课题进行数据分析),利用我们这学期所学知识建立恰当模型(ARIMA、GARCH 等),小组成员(不超过三人一组, 自由组队)中的一人安排 5 至 ...

    复习心烦,偶遇大作业,故摸鱼

    作业题目

    自由选取一组数据(可以是R 自带的数据集、或者其它来源,鼓励选取一些有趣的课题进行数据分析),利用我们这学期所学知识建立恰当模型(ARIMA、GARCH 等),小组成员(不超过三人一组, 自由组队)中的一人安排 5 至 10 分钟左右的PPT 课堂展示(教学周第16 周、无需展示代码,仅汇报初步结果即可,我给出点评后再完善后提交)。内容基本要求如下:
    • 说明数据来源背景, 需要分析解决的问题;
    • 包含完整的建模过程(拟合、估计、诊断检验、预测);
    • 结论或者建议(依赖于你拟解决的问题)。
    注意!!!数据量的选取需要适中, 留取部分原始数据作为评估模型预测能力使用。
    此部分对最后成绩的贡献比例为20%。

    Background of Data Set

    A multivariate yearly time series from 1960 to 1979 with variables (Included in dataset TSA)

    Test set and fit set

    Select the first 60 data for fitting , while the last 12 (one year) as the test data to test the final model we struct (cut the data frame by using R function window( ))

    Data preprocessing

    Observe the image directly and there is an obvious logarithmic upward trend. Carry out differential processing on dataset——wages to eliminate the upward trend,
    在这里插入图片描述
    test the unit root(DK test), the result are as followings:

    > adf.test(model)
    
    	Augmented Dickey-Fuller Test
    
    data:  model
    Dickey-Fuller = -3.6536, Lag order = 3, p-value = 0.03639
    alternative hypothesis: stationary
    

    We are not satisfied with the p p p value. In order to avoid excessive difference, we try to use Box-Cox transformation. From the results, it can be seen that the value of λ \lambda λ is 1. Although the value of p p p value for significance 0.01 is not good, it can only reduce the significance requirement ( α \alpha α = 0.05), rather than searching complex alternatives.
    在这里插入图片描述
    We finally preprocess the dataset by difference : model = diff (log(waves)). From the figure below, we can see that the peaks and troughs are distributed regularly. Combined with the actual background and common sense of life of the data set, workers’ salary will increase at the end of the year. After all, everyone wants to run. Workers’ salary should be fluctuated in an annual cycle. We can see from the image of data after difference. The salary at the end of each year is relatively high, but there are three hidden dangers:
    ① the periodicity is not very obvious. Although it has practical significance, for the size of sample is so small.
    ② Although it is expected that the salary at the end of middle age should be high, the data in 1986 conflicts with this statement, which is likely to be related to the social background at that time, and this information is unknown to us. The prediction work we do later may not be so good
    ③ It can also be seen from the figure that the growth rate of wages has slowed down (because we have just done differential processing, in fact, wages have been rising all the time)
    在这里插入图片描述

    Choose Model

    Calculate the ACF and PACF of the model. In fact, the diagrams of ACF and PACF can’t put forward feasible suggestions for us to choose p p p and q q q, hence it’s not necessary to calculate eacf. The results given must not be so good. We try to make a seasonal difference to the model.
    在这里插入图片描述
    After seasonal difference, it is also difficult to select appropriate models according to ACF and PACF images,
    在这里插入图片描述
    This model is a little troublesome, for model A R I M A ( p , d , q ) × ( P , D , Q ) s \rm ARIMA(p,d,q)\times(P,D,Q)_s ARIMA(p,d,q)×(P,D,Q)s,now we juat know the value of d d d and D D D A R I M A ( p , 1 , q ) × ( P , 1 , Q ) 11 ARIMA(p,1,q)\times(P,1,Q)_{11} ARIMA(p,1,q)×(P,1,Q)11, Next we select the optimal subset based on the BIC criterion. The following two figures respectively select the optimal subset from fit_set and model2.
    在这里插入图片描述
    We consider the models: A R I M A ( 0 , 1 , 0 ) × ( 0 , 1 , 1 ) 12 ARIMA(0,1,0)\times(0,1,1)_{12} ARIMA(0,1,0)×(0,1,1)12, A R I M A ( 0 , 1 , 0 ) × ( 1 , 1 , 1 ) 12 ARIMA(0,1,0)\times(1,1,1)_{12} ARIMA(0,1,0)×(1,1,1)12, A R I M A ( 0 , 1 , 1 ) × ( 0 , 1 , 0 ) 12 ARIMA(0,1,1)\times(0,1,0)_{12} ARIMA(0,1,1)×(0,1,0)12,
    A R I M A ( 0 , 1 , 0 ) × ( 1 , 1 , 0 ) 12 ARIMA(0,1,0)\times(1,1,0)_{12} ARIMA(0,1,0)×(1,1,0)12, A R I M A ( 1 , 1 , 0 ) × ( 0 , 1 , 0 ) 12 ARIMA(1,1,0)\times(0,1,0)_{12} ARIMA(1,1,0)×(0,1,0)12, A R I M A ( 1 , 1 , 0 ) × ( 1 , 1 , 1 ) 12 ARIMA(1,1,0)\times(1,1,1)_{12} ARIMA(1,1,0)×(1,1,1)12

    Parameter estimation

    a1=arima(fit_set,order=c(0,1,0),seasonal = list(order=c(0,1,1),period = 12))
    a2=arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,1),period = 12))
    a3=arima(fit_set,order=c(0,1,1),seasonal = list(order=c(0,1,0),period = 12))
    a4=arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,0),period = 12))
    a5=arima(fit_set,order=c(1,1,0),seasonal = list(order=c(0,1,0),period = 12))
    a6=arima(fit_set,order=c(1,1,0),seasonal = list(order=c(1,1,1),period = 12))
    

    Model diagnosis

    Diagnose the six models by R function tsdiag( ). During the diagnosing, it was found that the Ljung box test of model a1 was not good, one point was below the critical line and one point was very close to the critical line; The residual autocorrelation diagram of model a3 and model a5 is also not good, there’re some multiple outliers, so we exclude model a1, model a3 and model a5 first. There are still model a2, model a4 and model a6 models left. From the results of parameter estimation and model diagnosis, these models have no obvious advantages or disadvantages. However, based on the modeling principle of “simple but not simpler”, we exclude model a6 and leave model a2 and model a4 models for final prediction. The following are the results of testing a2 and a4 model.
    在这里插入图片描述
    在这里插入图片描述

    Model prediction

    a2=stats::arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,1),period = 12))
    a2_fore=forecast(a2,h=11,level = c(95))
    plot(a2_fore)
    lines(a2_fore$fitted,col="black")
    lines(wages,col="red")
    

    在这里插入图片描述

    a4=stats::arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,0),period = 12))
    a4_fore=forecast(a4,h=11,level = c(95))
    plot(a4_fore)
    lines(a4_fore$fitted,col="black")
    lines(wages,col="red")
    

    在这里插入图片描述

    Conclusion

    We can see the prediction is not successful, for each predictive value is a little higher than the true value. It may be caused by the logarithmic upward trend of the data. Differential processing is suitable for linear trends, but not for logarithmic trend. A possible solution is that we can add a approxiate penalization on t, in order to slow down the upward trend.

    Code

    library(forecast)
    library(TSA)
    library(xts)
    library(tseries)
    data(wages)
    help(wages)
    plot(wages)
    # 分组
    fit_set=window(wages,start=c(1981,7),end=c(1986,6))
    test_set=window(wages,start=c(1986,7),end=c(1987,6))
    # # 确定趋势
    # vec_wage=as.vector(wages)
    # newwage=data.frame(log_wage=log(vec_wage),day=c(1:length(vec_wage)))
    # fit=lm(newwage$log_wage~newwage$day,data=newwage)
    # summary(fit)
    # rm(fit)
    # 一次差分
    model=diff(fit_set)
    plot(model,ylab="model")
    abline(h=0)
    abline(v=c(1982:1986),lty=2)
    adf.test(model)
    # 尝试Box-Cox变换
    b=BoxCox.lambda(model)
    lambda=which.max(b)
    b=BoxCox.ar(model-min(model)*1.1)
    # 模型识别
    par(mfrow=c(1,2))
    acf(as.vector(model),lag.max=60,main="model")
    pacf(as.vector(model),lag.max = 60,main="model")
    # 考虑季节模型
    model2=diff(model,lag = 12)
    acf(as.vector(model2),lag.max=60,main="model2")
    pacf(as.vector(model2),lag.max = 60,main="model2")
    par(mfrow=c(1,1))
    # 模型识别
    par(mfrow=c(1,2))
    dis=armasubsets(y=fit_set,nar=13,nma=13,y.name="test",ar.method="ols")
    plot(dis)
    dis=armasubsets(y=model2,nar=4,nma=4,y.name="diff(1)",ar.method="ols")
    plot(dis)
    # 模型拟合
    a1=arima(fit_set,order=c(0,1,0),seasonal = list(order=c(0,1,1),period = 12))
    a2=arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,1),period = 12))
    a3=arima(fit_set,order=c(0,1,1),seasonal = list(order=c(0,1,0),period = 12))
    a4=arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,0),period = 12))
    a5=arima(fit_set,order=c(1,1,0),seasonal = list(order=c(0,1,0),period = 12))
    a6=arima(fit_set,order=c(1,1,0),seasonal = list(order=c(1,1,1),period = 12))
    # 模型诊断
    tsdiag(a1)
    tsdiag(a2)
    tsdiag(a3)
    tsdiag(a4)
    tsdiag(a5)
    tsdiag(a6)
    # 模型预测
    a2=stats::arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,1),period = 12))
    a2_fore=forecast(a2,h=11,level = c(95))
    plot(a2_fore)
    lines(a2_fore$fitted,col="black")
    lines(wages,col="red")
    
    a4=stats::arima(fit_set,order=c(0,1,0),seasonal = list(order=c(1,1,0),period = 12))
    a4_fore=forecast(a4,h=11,level = c(95))
    plot(a4_fore)
    lines(a4_fore$fitted,col="black")
    lines(wages,col="red")
    
    展开全文
  • 基于R语言时间序列分析所有指令[2021]

    千次阅读 多人点赞 2021-02-13 11:43:54
    1 安装包指令2 加载包指令3 help指令的使用4 读取不同格式数据4.1 读取csv格式的数据4.2 读取txt格式的数据4.3 读取xls和xlsx格式的数据4.4 参数使用5 ts生成时间序列的对象5.1 时间间隔为年的情况5.2 时间间隔为月...

        文章主要是总结一学期所学,基本覆盖了所有常见的指令,足够完成arima模型的数据选择到模型预测。
         时间序列应用广泛,不能仅仅局限于理论学习,代码实践更为重要。
    往期文章链接:
    基于 ARIMA-GARCH 模型人名币汇率分析与预测[论文完整][2020年]

    1 安装包指令

     install.packages("tseries")
     install.packages("forecast")
    

    2 加载包指令

    说明:载入一个包之后,就可以使用一系列新的函数和数据集了。包中往往提供了演示性的小型数据集和示例代码,能够让我们尝试这些新功能

     library(tseries)
     library(forecast)
     
    

    3 help指令的使用

    说明:使用函数help()可以查看其中任意函数或数据集的更多细节

    help(package="tseries")
    

    4 读取不同格式数据

    4.1 读取csv格式的数据

    data<-read.csv("路径",header=F,sep=",")
    

    4.2 读取txt格式的数据

    data<-read.table("路径",header=F,sep=",")
    

    4.3 读取xls和xlsx格式的数据

    可以统一将数据文件格式转成csv进行处理

    4.4 参数使用

    路径(注意反斜杠的方向)
    E:/大三上课程/时间序列.csv
    E:/大三上课程/时间序列.txt

    heaer
    这是一个逻辑值,T or F 反映文件的第一行是否包含变量名

    sep
    文件中的字段分隔符,例如对用制表符分割的文件使用sep="\t"

    5 ts生成时间序列的对象

    5.1 时间间隔为年的情况

    ts(data$V1,frequency = 1,start=c(1975))
    

    frequent=1表示以年为单位间隔 start说明了数据从1975年开始
    在这里插入图片描述

    5.2 时间间隔为月的情况

    ts(data$V1,frequency = 12,start=c(1975,1))
    

    frequent=12表示以月为单位间隔 start说明数据从1975年1月开始
    在这里插入图片描述

    5.3 时间间隔为日的情况

    ts(data$V1,frequency = 365,start=c(1975,1,1))
    

    frequent=365表示以日为单位间隔 start说明数据从1975年1月1日开始
    在这里插入图片描述

    5.4 参数使用

    data
    这个必须是一个矩阵,或者向量,再或者数据框frame

    补充 :data$v1 表示数据源的第一排 v2表示第二排 容易遗漏
    Frequency
    这个是时间观测频率数,也就是每个时间单位的数据数目
    Start
    时间序列开始值(是一个向量),允许第一个个时间单位出现数据缺失

    6 plot绘制时间序列的折线图

    plot基本用法

    plot(x=x轴数据,y=y轴数据,main="标题",sub="子标题",type="线型",xlab="x轴名称",ylab="y轴名称",xlim = c(x轴范围,x轴范围),ylim = c(y轴范围,y轴范围))
    

    实例演示

    ts_1<-ts(data$V1,frequency = 365,start=c(1975,1,1))
    plot(ts_1,main="ceshi",ylab="大小",xlab="时间",ylim=(-3:3))
    

    7 平稳性检验的三种方式

    7.1 图检法

            通过plot()绘制时间序列的时序图,如果序列是平稳的,那么序列应该是围绕一个均值上下随机波动。如果序列呈现递增递减则不是平稳时间序列。
            如下图时序图无明显趋势,在一个值上下波动,初步判断为平稳序列。
    在这里插入图片描述

    7.2 通过acf自相关系数图

    对于平稳时间序列,其自相关图一般随着阶数的递增,自相关系数会迅速衰减至0附近,而非平稳时间序列则可能存在先减后增或者周期性波动等变动。如下图所示,该时间序列随着阶数的递增迅速衰减至0附近,因此,可以判断该时间序列不是平稳时间序列。

    acf(ts_1)
    

    在这里插入图片描述

    7.3 通过adf函数p值判断(不需要主观判断)

    #使用前要加载包  ts_1是通过ts生成的时间序列(不清楚看前面)
    install.packages("tseries")
    library(tseries)
    adf.test(ts_1)
    

    在这里插入图片描述
    p值为0.02547小于显著性水平0.05接受备择假设,为平稳序列。若p值大于0.05则为非平稳性检验

    8 白噪声(纯随机)检验的方式

    #其中ts_1是通过ts做出的时序序列,for循环是为了多判断几期p值的判断。
    #一般我们判断前12期的p值,期数越大越不容易小于0.05for(i in 1:4print(Box.test(ts_1,lag=2*i))
    

    当p值都小于0.05,我们接受备择假设,该序列为非随机序列。
    当p值大于0.05,我们选择原假设,该序列为 纯随机序列。
    例如下,p值都小于0.05判断为非随机序列。如果为随机序列则数据之间无关联,只能进行差分或者放弃。

    在这里插入图片描述

    9 绘制自相关系数图和偏自相关系数图

    自相关系数图acf

    acf(ts_1)
    

    在这里插入图片描述
    偏自相关系数pacf

    pacf(ts_1)
    

    补充
    下面指令可以将直接做出时序图,acf,pacf

    tsdisplay(ts_1)
    

    在这里插入图片描述

    10模型选择及其定阶

    10.1通过acf图和pacf图

    • 若PACF序列满足在p步截尾,且ACF序列被负指数函数控制收敛到0,则Yn为AR§序列。

    • 若ACF序列满足在q步截尾,且PACF序列被负指数函数控制收敛到0,则Yn为MA(q)序列

    • 若ACF序列和PACF序列满足皆不截尾,但都被负指数函数控制收敛到0,则Yn为ARMA序列。

    当我们选择好模型可以通过如下指令进行构建模型

    #其中p表示AR的阶数 d表示差分的次数 q表示MA的阶数
    all<-arima(ts_1,order=c(p,d,q))
    all
    

    在这里插入图片描述

    10.2通过auto.arima()自动选择

    # 无特殊要求使用默认参数即可
    auto.arima(ts_1)
    

    如下函数为我们选择的是AR(1) ,aic aicc bic都是用来判断模型的好坏的准则,后面具体说明。
    在这里插入图片描述

    11AIC,BIC准则判断模型好坏

    AIC
    AIC = -2In(L) + 2k
    其中L指对应的最大似然函数,k指对应的模型的变量的个数。

    BIC
    BIC = -2In(L) + In(n)*k
    n指对应的数据数量,L和k同上所述。kln(n)惩罚项在维数过大且训练样本数据相对较少的情况下,可以有效避免出现维度灾难现象。

    AIC,BIC越小越好

    12 foreca进行模型预测

    # all表示模型 all<-arima(ts_1,order=c(1,0,0))
    # 60表示预测接下来的60个值
    library(forecast)
    plot(forecast(all,60))#绘制图形
    forecast(all,60) #显示所有预测值
    

    在这里插入图片描述

    13 差分命令

    #s1 表示差分后新的序列 
    #data$v1表示数据源的第一列,也就是要进行操作的数据
    #1表示差分一次 
    s1<-diff(data$v1,1)
    

    14 差分命令

    15 完整操作流程图

    在这里插入图片描述

    总结

    最后希望给文章点个赞,整理不易!!!
    最后希望给文章点个赞,整理不易!!!
    最后希望给文章点个赞,整理不易!!!
    如果有错误可以评论私信。
    如有需要完整论文及代码数据便于参考学习可评论、私信

    在这里插入图片描述

    展开全文
  • 时间序列分析及应用:R语言(原书第2版),市面上很好的一本书,pdf版,希望大家能受益
  • R语言时间序列分析教程学习,理论学习,代码学习。
  • 蔡瑞胸的多元时间序列分析及金融应用 R语言>>中文版,带书签
  • 王燕。时间序列分析 R语言 案例和习题,用于R语言的操作练习
  • 时间序列数据分析 R软件应用.pdf
  • 使用R语言进行时间序列分析

    千次阅读 2020-04-15 15:10:54
    时间序列是将统一统计值按照时间发生的先后顺序来进行排列,时间序列分析的主要目的是根据已有数据对未来进行预测。 一个稳定的时间序列中常常包含两个部分,那么就是:有规律的时间序列+噪声。所以,在以下的方法中...
  • 本书以易于理解的方式讲述了时间序列模型及其应用,内容包括:趋势,平稳时间序列模型,非平稳时间序列模型,模型识别、参数估计,模型诊断,预测、季节模型、时间序列回归
  • 时间序列分析及应用:R语言 原书第2版 》以易于理解的方式讲述了时间序列模型及其应用 主要内容包括:趋势 平稳时间序列模型 非平稳时间序列模型 模型识别 参数估计 模型诊断 预测 季节模型 时间序列回归模型 异...
  • 里面包含数据集,R语言代码,代码后面也写了注释,清晰易懂
  • 一个进行时序分析的实例,所使用的模型为乘法季节ARIMA模型
  • 主要使用R语言时间序列分析,详细介绍时间序列如何进行分析以及r语言涉及其中的代码。
  • R语言时间序列分析实例

    千次阅读 2020-12-22 12:50:40
    #加载数据x=read.table(file.choose())#生成时间序列对象xtimeseries#画时间序列图plot.ts(xtimeseries)#增加线性拟合曲线abline(lm(xtimeseries~time(xtimeseries)))1、分解时间序列分解一个时间序列意味着把它拆分...
  • 本文是我个人对于之前学习的R语言的一个复习,主要目的是便于理解和使用,文中并没有过多对于函数原理计算公式的介绍,主要是何处应用和怎样应用,如果对于具体的原理感兴趣,可以查阅书籍或自行网上搜索。
  • R语言时间序列练习

    2019-04-19 15:49:36
    R语言时间序列数据分析,常用的TTR包的导入、入门与作图。
  • 拿到的GIMMS NDVI数据是ENVI格式(不是NASA的原始数据),要求做时间序列分析,查看变化情况,那么该怎么做呢?
  • 黄天元,复旦大学博士在读,热爱数据科学与开源工具(R),致力于利用数据科学迅速积累行业经验优势和科学知识发现,涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等,著有《R语言数据...
  • 时间序列分析及应用:R语言

    热门讨论 2014-08-12 10:44:41
    时间序列分析及应用:R语言,(原书第2版)。内容:应用R语言实现时间序列分析的具体方法与步骤,很好的一本学习用书。
  • 复旦大学博士在读,热爱数据科学与开源工具(R),致力于利用数据科学迅速积累行业经验优势和科学知识发现,涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等,著有《R语言数据高效...
  • 在上一篇中,探讨了R语言时间序列分析常用步骤,如何比对AIC值判断最优模型?代码和解释如下: #WWWusage是datasets包自带的每分钟通过服务器连接到因特网的用户数的长度为100的时间序列数据 require(graphics) #...
  • 用以学习时间序列分析,书中有相应的R程序,以易于理解的方式讲述了时间序列模型及其应用,内容包括趋势、平稳时间序列模型、非平稳时间序列模型、模型识别、模型诊断、预测等等

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 68,273
精华内容 27,309
关键字:

r语言 时间序列分析