8时间序列回归模型——R实现

时间序列回归模型

1干预分析

1.1概念及模型

Box和Tiao引入的干预分析提供了对于干预影响时间序列的效果进行评估的一个框架，假设干预是可以通过时间序列的均值函数或者趋势而对过程施加影响，干预可以自然产生也可以人为施加的，如国家的宏观调控等。

其模型可以如下表示：

其中mt代表均值的变化，Nt是ARIMA过程。

1.2干预的分类

阶梯响应干预

脉冲响应干预

1.3干预的实例分析

1.3.1模型初探

对数化航空客运里程的干预模型的估计

> data(airmiles)

> acf(as.vector(diff(diff(window(log(airmiles),end=c(2001,8)),12))),lag.max=48)#用window得到在911事件以前的未爱干预的时间序列子集

对暂用的模型进行诊断

>fitmode<-arima(airmiles,order=c(0,1,1),seasonal=list(order=c(0,1,0)))

> tsdiag(fitmode)

从诊断图可以看出存在三个异常点，acf在12阶存在高度相关因此在季节中加入MA（1）系数。

1.3.2拟合带有干预信息的模型

函数：

arimax(x, order = c(0, 0, 0), seasonal = list(order = c(0, 0, 0), period = NA), xreg = NULL, include.mean = TRUE, transform.pars = TRUE, fixed = NULL, init = NULL, method = c("CSS-ML", "ML", "CSS"), n.cond, optim.control = list(), kappa = 1e+06, io = NULL, xtransf, transfer = NULL)

arimax函数扩展了arima函数，可以处理时间序列中干扰分析及异常值。假设干扰影响过程的均值，相对未受干扰的无价值函数的偏离用一些协变量的ARMA滤波器的输出这种来表示，偏差被称作传递函数。构造传递函数的协变量通过xtransf参数以矩阵或者data.frame 的形式代入arimax函数。

air.m1=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1),

period=12),xtransf=data.frame(I911=1*(seq(airmiles)==69),

I911=1*(seq(airmiles)==69)),

transfer=list(c(0,0),c(1,0)),xreg=data.frame(Dec96=1*(seq(airmiles)==12),

Jan97=1*(seq(airmiles)==13),Dec02=1*(seq(airmiles)==84)),method='ML')

> air.m1

Call:

arimax(x = log(airmiles), order = c(0, 1, 1), seasonal = list(order = c(0, 1,

1), period = 12), xreg = data.frame(Dec96 = 1 * (seq(airmiles) == 12), Jan97 = 1 *

(seq(airmiles) == 13), Dec02 = 1 * (seq(airmiles) == 84)), method = "ML",

xtransf = data.frame(I911 = 1 * (seq(airmiles) == 69), I911 = 1 * (seq(airmiles) ==

69)), transfer = list(c(0, 0), c(1, 0)))

Coefficients:

ma1 sma1 Dec96 Jan97 Dec02 I911-MA0 I911.1-AR1 I911.1-MA0

-0.3825 -0.6499 0.0989 -0.0690 0.0810 -0.0949 0.8139 -0.2715

s.e. 0.0926 0.1189 0.0228 0.0218 0.0202 0.0462 0.0978 0.0439

sigma^2 estimated as 0.0006721: log likelihood = 219.99, aic = -423.98

画图

plot(log(airmiles),ylab="log(airmiles)")

points(fitted(air.m1))

Nine11p=1*(seq(airmiles)==69)

plot(ts(Nine11p*(-0.0949)+

filter(Nine11p,filter=.8139,method='recursive',side=1)*(-0.2715),

frequency=12,start=1996),type='h',ylab='9/11 Effects')

abline(h=0)

从上图可以看出在2003年底后，911事件的影响效应才平息，航班客运量恢复了正常。2异常值

在时间序列中异常有两种，可加异常和新息异常，分别记AO和IO。

2.1异常值示例

2.1.1模拟数据

模拟一般的ARIMA（1，0，1），然后故意将第10个观测值变成异常值10.

> set.seed(12345)

> y=arima.sim(model=list(ar=0.8,ma=0.5),n.start=158,n=100)

> y

Time Series:

Start = 1

End = 100

Frequency = 1

[1] 0.49180881 -0.22323665 -0.99151270 -0.73387818 -0.67750094

-1.14472133 -2.14844671 -2.49530794