zoukankan      html  css  js  c++  java
  • R语言时间序列ARMAX建模

    建立ARMAX模型需要运用R的dse包,在R的dse包中The ARMA model representation is general, so that VAR, VARX,ARIMA, ARMAX, ARIMAX can all be considered to be special cases.

    数据集为天然气炉中的天然气(input)与产生的CO2(output),数据来源为王燕应用时间序列分析第三版附录表A1-24,首先队数据做简要的分析,做出时序图以及协方差图

    input<-ts(input,frequency = 1,start = c(1,1))
    
    output<-ts(output,frequency = 1,start = c(1,1))
    
    par(mfrow=c(1,2))
    
    plot.ts(input)
    
    plot.ts(output)
    

      

    做两个变量间的协方差函数需要使用ccf函数

    ccf(output,input,lag.max = NULL,type = c("correlation"),plot = TRUE)
    

      

    可以看到两变量存在逻辑上的因果关系,变量延迟的阶数在3~5左右,因此可以通过dse包进行建模,需要使用的函数有ARMA,TSdata,estMaxLik,forecast

    首先使用ARMA函数定义一个模型表达式,注意ARMA只能定义一个模型结构,而不能估计其参数,其定义的结构为

    A(L)y(t) = B(L)w(t) + C(L)u(t) + TREND(t)

    where

    A (axpxp) is the auto-regressive polynomial array.
    B (bxpxp) is the moving-average polynomial array.
    C (cxpxm) is the input polynomial array. C should be NULL if there is no input
    y is the p dimensional output data.
    u is the m dimensional control (input) data.

    TREND is a matrix the same dimension as y, a p-vector (which gets replicated for each time
    period), or NULL.

    需要注意定义模型中参数与变量的区别,以及模型设置的初始值设置

    By default, elements in parameter arrays are treated as constants if they are exactly 1.0 or 0.0, and
    as parameters otherwise. A value of 1.001 would be treated as a parameter, and this is the easiest
    way to initialize an element which is not to be treated as a constant of value 1.0. Any array elements
    can be fixed to constants by specifying the list constants. Arrays which are not specified in the list
    will be treated in the default way. An alternative for fixing constants is the function fixConstants.

    因此我们可以如下定义模型

    MA <- array(c(1,-1.5,0.5),c(3,1,1))
    
    C <- array(c(0,0,0,-0.5,-0.5,-0.5),c(5,1,1))
    
    AR <- array(c(1,-0.5),c(2,1,1))
    
    TR <-array(c(50),c(1,1,1))
    
    ARMAX <- ARMA(A=AR, B=MA, C=C ,TREND<-TR)
    

      

    定义好模型后,还需要将数据转化为需要的格式,这里要使用TSdata函数,把input与output数据存储在一起

    随后就可以使用estMaxLik函数对模型进行求解了,先说另一个求解函数estVARXls,顾名思义,采用的求解方法为ls,该函数的好处是可以自动定阶,如:

    armaxdata<-TSdata(input=input,output=output)
    
    model <- estVARXls(armaxdata)
    

      

    > model
    neg. log likelihood= 821.601

    A(L) =
    1-1.705922L1+0.6411281L2+0.2290077L3-0.128453L4-0.0686206L5+0.03279703L6

    B(L) =
    1

    C(L) =
    -0.08792481+0.2803458L1-0.2801545L2-0.3750078L3-0.07263863L4+0.5139702L5

    如果想要用前面设置的ARMA形式来建模的话则需要使用estMaxLik,该函数采用Maximum likelihood estimation,注意每次计算的结果可能与ARMA设置的初始值有关,如果想要得到更好的结果,需要对该函数参数进行讨论。

    model <- estMaxLik(armaxdata, ARMAX, algorithm="optim")
    

      

    > model
    neg. log likelihood= 231.8164
    TREND= 53.493
    A(L) =
    1+0.001700993L1

    B(L) =
    1+1.486967L1+0.7612954L2

    C(L) =
    0-1.002061L3-1.601826L4

    最后还要使用forecast函数对模型进行预测,注意如果input和output一样长,则不会产生新的未来预测值,但是会产生模型拟合值pred,如果想要对未来进行预测,需要讨论conditioning.inputs.forecasts(A time series matrix or list of time series matrices to append to input variables for the forecast periods.)

    pr <- forecast(model, conditioning.inputs=inputData(armaxdata))
    

      

    提取出pred项画出预测图

    pred<-ts(as.vector(pr$pred),frequency = 1,start = c(1,1))
    plot(pred,col="red",type="p")
    par(new=TRUE)
    plot(output)
    

      

     笔者也是ARMAX初学者,若有错误请大家指出

    博主原创内容,禁止一切形式转载
  • 相关阅读:
    使用jsonp跨域调用百度js实现搜索框智能提示(转)
    jsonp 跨域
    Aixs2发布webservice服务
    java web service 上传下载文件
    java 网页 保存上传文件
    flash、js 函数 互相调用
    java web工程启动socket服务
    mysql 在Windows下自动备份
    Servlet中几个常用方法的推衍
    Tomcat常用设置 <持续更新>
  • 原文地址:https://www.cnblogs.com/statruidong/p/6910893.html
Copyright © 2011-2022 走看看