时间序列分析
This presentation is the property of its rightful owner.
Sponsored Links
1 / 141

时间序列分析 PowerPoint PPT Presentation


  • 301 Views
  • Uploaded on
  • Presentation posted in: General

时间序列分析. 考虑 …. 何为科学 ? 何为统计 ? 统计是科学吗 ? 什么是数学 ? 统计是数学吗 ? 能够证明某些模型或理论是正确的吗 ?. 统计. 归纳 : 从部分到总体,从特殊到一般的推理 . 演绎 : 基于假定,公理的严格的证明链. 时间序列的例子。从图形你会想到 : 数据会是什么样子的 , 可能的模式及它们的意义 , 可能的预测,如何点图 ,…….

Download Presentation

时间序列分析

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


5555123

时间序列分析


5555123

考虑…

  • 何为科学?

  • 何为统计?

  • 统计是科学吗?

  • 什么是数学?

  • 统计是数学吗?

  • 能够证明某些模型或理论是正确的吗?

Stat 153, Xizhi Wu


5555123

统计

  • 归纳: 从部分到总体,从特殊到一般的推理.

  • 演绎: 基于假定,公理的严格的证明链.

Stat 153, Xizhi Wu


5555123

时间序列的例子。从图形你会想到: 数据会是什么样子的, 可能的模式及它们的意义, 可能的预测,如何点图,……

Stat 153, Xizhi Wu


Examples economic and financial time series beverage wheat price annual index series from 1500 1864

Examples: Economic and financial time series Beverage wheat price annual index series from 1500-1864

Stat 153, Xizhi Wu

T01.R


Examples economic and financial time series beverage wheat price annual index series from 1810 1864

Examples: Economic and financial time series Beverage wheat price annual index series from 1810-1864

Stat 153, Xizhi Wu

T01.R


5555123

Examples: Physical time seriesAverage air temperature (deg C) at Recife, Brazil, in successive months from 1953-1962

Stat 153, Xizhi Wu

T01.R


5555123

Examples: Marketing time series Sales of an industrial beater in successive months from January 1965 to November 1971.

Stat 153, Xizhi Wu

T01.R


Examples f i nance the percentage of beijing residents long term deposit over the total balance

Examples: Finance The Percentage of Beijing Residents' Long-term Deposit Over the Total Balance

Stat 153, Xizhi Wu

T01.R


Examples pr o cess control data

Examples: Process control data

Stat 153, Xizhi Wu

T01.R


Examples binary process

Examples: Binary process

Stat 153, Xizhi Wu

T01.R


Examples point process

Examples: Point process

Stat 153, Xizhi Wu

T01.R


5555123

你有时间序列的例子吗?

Stat 153, Xizhi Wu


5555123

术语

  • 连续时间序列

  • 离散时间序列

  • 通常是等间距的:

  • 抽样的序列

  • 即时的或整合的

  • 通常不是独立的

  • 确定的或随机的(精确预测是不可能的)

Stat 153, Xizhi Wu


Passengers in china s airports

看看这个时间序列Passengers in China’s Airports

Stat 153, Xizhi Wu


Passengers in china s airports1

看看这个时间序列Passengers in China’s Airports

Stat 153, Xizhi Wu


5555123

横截面数据时间序列数据

  • 人们对统计数据往往可以根据其特点从两个方面来切入,以简化分析过程。一个是研究所谓横截面(cross section)数据,也就是对大体上同时,或者和时间无关的不同对象的观测值组成的数据。

  • 另一个称为时间序列(time series),也就是由对象在不同时间的观测值形成的数据。

  • 前面讨论的模型多是和横截面数据有关。这里将讨论时间序列的分析。我们将不讨论更加复杂的包含这两方面的数据。


5555123

时间序列和回归

  • 时间序列分析也是一种回归。

  • 回归分析的目的是建立因变量和自变量之间关系的模型;并且可以用自变量来对因变量进行预测。通常线性回归分析因变量的观测值假定是互相独立并且有同样分布。

  • 而时间序列的最大特点是观测值并不独立。时间序列的一个目的是用变量过去的观测值来预测同一变量的未来值。也就是说,时间序列的因变量为变量未来的可能值,而用来预测的自变量中就包含该变量的一系列历史观测值。

  • 当然时间序列的自变量也可能包含随着时间度量的独立变量。


15 1 tax txt tax sav

例15.1 (数据:Tax.txt,Tax.sav)

某地从1995年1月到2005年7月的税收(单位:万元)。该数据有按照时间顺序的按月记录,共127个观测值。图15.1就是由该数据得到的一个时间序列图。


15 1 tax txt tax sav1

例15.1 (数据:Tax.txt,Tax.sav)

从这个点图可以看出。总的趋势是增长的,但增长并不是单调上升的;有涨有落。大体上看,这种升降不是杂乱无章的,和季节或月份的周期有关系。当然,除了增长的趋势和季节影响之外,还有些无规律的随机因素的作用。这个只有一种随着时间变化的变量(税收)的序列一般称为纯粹时间序列(pure time series)。下面将通过该例子对纯粹时间序列进行介绍。


5555123

时间序列的组成部分

  • 从该例可以看出,该时间序列可以有三部分组成:趋势(trend)、季节(seasonal)成分和无法用趋势和季节模式解释的随机干扰(disturbance)。

  • 例中数据的税收就就可以用这三个成分叠加而成的模型来描述。

  • 一般的时间序列还可能有循环或波动(Cyclic, or fluctuations)成分;循环模式和有规律的季节模式不同,周期长短不一定固定。比如经济危机周期,金融危机周期等等。


5555123

时间序列的组成部分

  • 一个时间序列可能有趋势、季节、循环这三个成分中的某些或全部再加上随机成分。因此,

  • 如果要想对一个时间序列本身进行较深入的研究,把序列的这些成分分解出来、或者把它们过虑掉则会有很大的帮助。

  • 如果要进行预测,则最好把模型中的与这些成分有关的参数估计出来。

  • 就例中的时间序列的分解,通过SPSS软件,可以很轻而易举地得到该序列的趋势、季节和误差成分。


5555123

去掉季节成分,只有趋势和误差成分的例15.1的时间序列。


5555123

例15.1的时间序列分解出来的纯趋势成分和纯季节成分两条曲线


5555123

例15.1的时间序列分解出来的纯趋势成分和纯误差成分两条曲线


5555123

指数平滑

  • 如果我们不仅仅满足于分解现有的时间序列,而且想要对未来进行预测,就需要建立模型。首先,这里介绍比较简单的指数平滑(exponential smoothing)。

  • 指数平滑只能用于纯粹时间序列的情况,而不能用于含有独立变量时间序列的因果关系的研究。

  • 指数平滑的原理为:当利用过去观测值的加权平均来预测未来的观测值时(这个过程称为平滑),离得越近的观测值要给以更多的权。

  • 而“指数”意味着:按照已有观测值“老”的程度,其上的权数按指数速度递减。


5555123

指数平滑

  • 以简单的没有趋势和没有季节成分的纯粹时间序列为例,指数平滑在数学上这实际上是一个几何级数。这时,如果用Yt表示在t时间的平滑后的数据(或预测值),而用X1, X2, …, Xt表示原始的时间序列。那么指数平滑模型为

或者,等价地,

这里的系数为几何级数。因此称之为“几何平滑”比使人不解的“指数平滑”似乎更有道理。


5555123

指数平滑

  • 自然,这种在简单情况下导出的公式(如上面的公式)无法应对具有各种成分的复杂情况。

  • 后面将给出各种实用的指数平滑模型的公式。

  • 根据数据,可以得到这些模型参数的估计以及对未来的预测。在和我们例子有关的指数平滑模型中,需要估计12个季节指标和三个参数(包含前面公式权重中的a,和趋势有关的g,以及和季节指标有关的d)。

  • 在简单的选项之后,SPSS通过指数平滑产生了对2005年6月后一年的预测。下图为原始的时间序列和预测的时间序列(光滑后的)。下面为误差。


5555123

例15.1时间序列数据的指数平滑和对未来的预测


5555123

x=scan("d:/booktj1/data/tax.txt")

tax=ts(x, frequency = 12, start = c(1995, 1))

ts.plot(tax,ylab="Tax")

#plot(x1,ylab="Sales")

a=stl(tax, "period") #分解

a$time.series #分解结果(三列)

ts.plot(a$time.series[,1:3])

b=HoltWinters(tax,beta=0) #Holt-Winters滤波指数平滑

predict(b,n.ahead=12) #对未来12个月预测

pacf(tax); acf(tax)

w=arima(tax, c(0, 1, 1),seasonal = list(order=c(1,2 ,1), period=12))

predict(w, n.ahead = 12)

w$residuals#残差

acf(w$resi)

pacf(w$resi)

w$coef#估计的模型系数

w$aic #aic值


Box jenkins arima

Box-Jenkins 方法:ARIMA模型

  • 如果要对比较复杂的纯粹时间序列进行细致的分析,指数平滑往往是无法满足要求的.

  • 而若想对有独立变量的时间序列进行预测,指数平滑更是无能为力。

  • 于是需要更加强有力的模型。这就是下面要介绍的Box-Jenkins ARIMA模型。

  • 数学上,指数平滑仅仅是ARIMA模型的特例.


Arima ar

ARIMA模型 :AR模型

  • 比指数平滑要有用和精细得多的模型是Box-Jenkins引入的ARIMA模型。或称为整合自回归移动平均模型(ARIMA 为Autoregressive Integrated Moving Average一些关键字母的缩写)。该模型的基础是自回归和移动平均模型或ARMA(Autoregressive and Moving Average) 模型。

  • 它由两个特殊模型发展而成,一个特例是自回归模型或AR (Autoregressive) 模型。假定时间序列用X1, X2, …, Xt表示,则一个纯粹的AR (p)模型意味着变量的一个观测值由其以前的p个观测值的线性组合加上随机误差项at(该误差为独立无关的)而得:

这看上去象自己对自己回归一样,所以称为自回归模型;它牵涉到过去p个观测值(相关的观测值间隔最多为p个.


Arima ma

ARIMA模型 :MA模型

  • ARMA模型的另一个特例为移动平均模型或MA (Moving Average) 模型,一个纯粹的MA (q)模型意味着变量的一个观测值由目前的和先前的q个随机误差的线性的组合:

由于右边系数的和不为1(q甚至不一定是正数),因此叫做“移动平均”不如叫做“移动线性组合”更确切;虽然行家已经习惯于叫“平均”了,但初学者还是因此可能和初等平滑方法中的什么“三点平均”之类的术语混淆。


Arima arma

ARIMA模型 :ARMA模型

  • 显然ARMA(p,q)模型应该为AR (p)模型和MA(q)模型的组合了:

显然ARMA(p,0)模型就是AR (p)模型,而ARMA(0,q)模型就是MA(q)模型。这个一般模型有p+q个参数要估计,看起来很繁琐,但利用计算机软件则是常规运算;并不复杂。


Arima

ARIMA模型:平稳性和可逆性

  • 但是要想ARMA(p,q)模型有意义则要求时间序列满足平稳性(stationarity)和可逆性(invertibility)的条件,

  • 这意味着序列均值不随着时间增加或减少,序列的方差不随时间变化,另外序列本身相关的模式不改变等。

  • 一个实际的时间序列是否满足这些条件是无法在数学上验证的,

  • 这没有关系,但可以从下面要介绍的时间序列的自相关函数和偏相关函数图中可以识别出来。

  • 一般人们所关注的的有趋势和季节/循环成分的时间序列都不是平稳的。这时就需要对时间序列进行差分(difference)来消除这些使序列不平稳的成分,而使其变成平稳的时间序列,并估计ARMA模型,估计之后再转变该模型,使之适应于差分之前的序列(这个过程和差分相反,所以称为整合的(integrated)ARMA模型),得到的模型于是称为ARIMA模型。


Arima1

ARIMA模型:差分

  • 差分是什么意思呢?差分可以是每一个观测值减去其前面的一个观测值,即Xt-Xt-1。这样,如果时间序列有一个斜率不变的趋势,经过这样的差分之后,该趋势就会被消除了。

  • 当然差分也可以是每一个观测值减去其前面任意间隔的一个观测值;比如存在周期固定为s的季节成分,

  • 那么相隔s的差分 为Xt-Xt-s就可以把这种以s为周期的季节成分消除。

  • 对于复杂情况,可能要进行多次差分,才能够使得变换后的时间序列平稳。


5555123

ARMA模型的识别和估计

  • 上面引进了一些必要的术语和概念。下面就如何识别模型进行说明。

  • 要想拟合ARIMA模型,必须先把它利用差分变成ARMA(p,q)模型,并确定是否平稳,然后确定参数p,q。

  • 现在利用一个例子来说明如何识别一个AR(p)模型和参数p。

  • 由此MA(q)及ARMA(p,q)模型模型可用类似的方法来识别。


5555123

ARMA模型的识别和估计

  • 根据ARMA(p,q)模型的定义,它的参数p,q和自相关函数(acf,autocorrelations function)及偏自相关函数(pacf,partial autocorrelations function)有关。

  • 自相关函数描述观测值和前面的观测值的相关系数;

  • 而偏自相关函数为在给定中间观测值的条件下观测值和前面某间隔的观测值的相关系数。

  • 这里当然不打算讨论这两个概念的细节。引进这两个概念主要是为了能够了解如何通过研究关于这两个函数的acf和pacf图来识别模型。


Ar2 sav

例:数据AR2.sav为了直观地理解上面的概念,下面利用一个数据例子来描述。


Ar2 sav acf pacf

例:数据AR2.sav;拖尾和截尾先来看该时间序列的acf(左)和pacf图(右)

左边的acf条形图是衰减的指数型的波动;这种图形称为拖尾。而右边的pacf条形图是在第二个条(p=2)之后就很小,而且没有什么模式;这种图形称为在在p=2后截尾。这说明该数据满足是平稳的AR(2)模型。


5555123

拖尾和截尾

  • 所谓拖尾图形模式也可能不是以指数形式,而是以正负相间的正弦形式衰减。类似地,如果acf图形是在第q=k个条后截尾,而pacf图形为拖尾,则数据满足MA(q)模型。如果两个图形都拖尾则可能满足ARMA(p,q)模型。具体判别法总结在下面表中(并不一定严格!):


Acf pacf

acf和pacf图

  • 如acf和pacf图中至少一个不是以指数形式或正弦形式衰减,那么说明该序列不是平稳序列,必须进行差分变换来得到一个可以估计参数的满足ARMA(p,q)模型的序列

  • 如一个时间序列的acf和pacf图没有任何模式,而且数值很小,那么这个序列可能就是一些互相独立的无关的随机变量。一个很好拟合的时间序列模型的残差就应该有这样的acf和pacf图。


Ar 2 ma 2 arma 2 2 acf pacf 0 p q

AR(2)、MA(2)和ARMA(2,2)模型所对应的acf和pacf图。注意,图中有些条是从0开始的(不算在p或q内)。


Ar2 sav1

例:数据AR2.sav

  • 根据acf和pacf图的形态,不用进行任何差分就可以直接用AR(2)模型拟合。利用SPSS软件,选择AR(2)模型(等价地ARIMA(2,0,0)(0,0,0)模型),得到参数估计为

也就是说该AR(2)模型为


Ar2 sav2

例:数据AR2.sav

  • 例15.2的原始序列和由模型AR(2)得到的拟合值及对未来10个观测的预测图(虚线)


Ar2 sav3

例:数据AR2.sav

  • 下面再看剩下的残差序列是否还有什么模式。这还可以由残差的pacf(左)和acf(右)图来判断。可以看出,它们没有什么模式;这说明拟合比较成功。


Ar2 sav4

例:数据AR2.sav

  • 下图为残差对拟合值的散点图。看不出任何模式。说明残差的确是独立的和随机的。


Arima p d q p d q s

ARIMA (p,d,q)(P,D,Q)s模型

  • 在对含有季节和趋势/循环等成分的时间序列进行ARIMA模型的拟合研究和预测时,就不象对纯粹的满足可解条件的ARMA模型那么简单了。

  • 一般的ARIMA模型有多个参数,没有季节成分的可以记为ARIMA(p,d,q),如果没有必要利用差分来消除趋势或循环成分时,差分阶数d=0,模型为ARIMA(p,0,q),即ARMA(p, q)。

  • 在有已知的固定周期s时,模型多了4个参数,可记为ARIMA(p,d,q)(P,D,Q)s。这里增加的除了周期s已知之外,还有描述季节本身的ARIMA(P,D,Q)的模型识别问题。因此,实际建模要复杂得多。需要经过反复比较。


Arima tax sav

用ARIMA模型拟合tax.sav

  • 先前对例15.1(数据tax.txt或tax.sav)进行了分解,并且用指数平滑做了预测。知道其中有季节和趋势成分。

  • 下面试图对其进行ARIMA模型拟合。先试图对该序列做acf和pacf条形图。其中acf图显然不是拖尾(不是以指数速率递减),因此说明需要进行差分。


Arima tax sav1

用ARIMA模型拟合例tax.sav

  • 关于于参数,不要选得过大;每次拟合之后要检查残差的acf和pacf图,看是否为无关随机序列。

  • 在SPSS软件中还有类似于回归系数的检验以及其他一些判别标准的计算机输出可做参考(这里不细说)。

  • 经过几次对比之后,对于例16.1数据我们最后选中了ARIMA(0,1,1)(1,2,1)12模型来拟合。拟合的结果和对2005年6月之后12个月的预测在下图中


5555123

例tax.sav的原始序列和由模型得到的拟合值及对未来12个月的预测图。


Tax sav

例:数据tax.sav

  • 为了核对,当然要画出残差的acf和pacf的条形图来看是否还有什么非随机的因素存在。下图为这两个点图,看来我们的模型选择还是适当的。


Tax sav1

例:数据tax.sav

  • 例15.1数据拟合ARIMA(0,1,1)(1,2,1)12模型后残差序列的Ljung-Box检验的p值


5555123

新SPSS


Arima2

用ARIMA模型拟合带有独立变量的时间序列

  • 例:数据tsadds2.sav是一个销售时间序列,以每周七天为一个季节周期,除了销售额序列sales之外,还有一个广告花费的独立变量adds。先不理睬这个独立变量,把该序列当成纯粹时间序列来用ARIMA模型拟合。右图为该序列的点图。


Tsadds2 sav

数据tsadds2.sav

  • 再首先点出其acf和pacf条形图

acf图显然不是拖尾模式,因此,必须进行差分以消除季节影响。试验多次之后,看上去ARIMA(2,1,2)( 0,1,1)7的结果还可以接受。残差的pacf和acf条形图在下一页图中


Arima3

用ARIMA模型拟合带有独立变量的时间序列

  • 继续改进我们的模型,再把独立变量广告支出加入模型,最后得到的带有独立变量adds的ARIMA(2,1,2)( 0,1,1)7模型。拟合后的残差图在下图中。


Arima4

用ARIMA模型拟合带有独立变量的时间序列

  • 从各种角度来看拟合带独立变量平方的ARIMA(2,1,2)( 0,1,1)7模型给出更好的结果。

  • 虽然从上面的acf和pacf图看不出(一般也不应该看出)独立变量对序列的自相关性的影响,但是根据另外的一些判别准则,独立变量的影响是显著的,而且加入独立变量使得模型更加有效。


Arima5

用ARIMA模型拟合带有独立变量的时间序列

  • 要注意,一些独立变量的效果也可能是满足某些时间序列模型的,也可能会和季节、趋势等效应混杂起来不易分辩。这时,模型选择可能就比较困难。也可能不同模型会有类似的效果。

  • 一个时间序列在各种相关的因素影响下的模型选择并不是一件简单明了的事情。实际上没有任何统计模型是绝对正确的,它们的区别在于,在某种意义上,一些模型的某些性质可能要优于另外一些。


5555123

新SPSS的时间序列实现

  • 特点:

    • 在ARIMA中自动选择用什么参数

    • 在指数平滑和ARIMA中自动选择模型(包括参数)

  • 下面是两个例子

    • TAX

    • AIRPORT


Tsadds2 sav1

tsadds2.sav


Airport sav

Airport.sav


5555123

GET

FILE='C:\XZWU\HEPbook\data\Airport.sav'.

DATASET NAME DataSet1 WINDOW=FRONT.

PREDICT THRU END.

* Time Series Modeler.

TSMODEL

/MODELSUMMARY PRINT=[ MODELFIT]

/MODELSTATISTICS DISPLAY=YES MODELFIT=[ SRSQUARE]

/SERIESPLOT OBSERVED FORECAST

/OUTPUTFILTER DISPLAY=ALLMODELS

/AUXILIARY CILEVEL=95 MAXACFLAGS=24

/MISSING USERMISSING=EXCLUDE

/MODEL DEPENDENT=北京 成都 福州 桂林 杭州 济南 老河口 南京 上海 武汉

盐城 银川

PREFIX='Model'

/EXPERTMODELER TYPE=[ARIMA EXSMOOTH] TRYSEASONAL=YES

/AUTOOUTLIER DETECT=OFF .


5555123

PREDICT THRU YEAR 2004 .

* Time Series Modeler.

TSMODEL

/MODELSUMMARY PRINT=[ MODELFIT]

/MODELSTATISTICS DISPLAY=YES MODELFIT=[ SRSQUARE]

/SERIESPLOT OBSERVED FORECAST FIT

/OUTPUTFILTER DISPLAY=ALLMODELS

/SAVE PREDICTED(Predicted) LCL(LCL) UCL(UCL)

/AUXILIARY CILEVEL=95 MAXACFLAGS=24

/MISSING USERMISSING=EXCLUDE

/MODEL DEPENDENT=北京 成都 福州 桂林 杭州 济南 老河口 南京 上海 武汉

盐城 银川

PREFIX='Model'

/EXPERTMODELER TYPE=[ARIMA EXSMOOTH] TRYSEASONAL=YES

/AUTOOUTLIER DETECT=OFF .


5555123

状态空间的计算例子

吴喜之


Model shumway notation transition equation and observation equation

Model (Shumway notation)Transition equation and observation equation

  • xt : p×1 state vector N(m, S)

  • yt : q×1 observation/measurement vector

  • wt : p×1 system noise iid N(0,Q)

  • vt : q×1 observation noise iid N(0,R)

  • ut : r×1 fixed input

  • At : q×p observation/measurement matrix

  • F: p×p ; G: q×r, U: p×r, Q: p×p, R: q×q


5555123

在骨移植后的白细胞数的纵向数据

#A Biometrical Example

x=scan("c:/xzwu/2009/berkeley/ts/shumway/mydata/Jones.dat",na="0")

y=matrix(x,91,3);

y=as.data.frame(y)

names(y)=c("log(white blood count)","log(platelet)","hematocrit(HCT)")

par(mfrow=c(3,1))

plot(ts(y[,1]),type="o",ylab="LWC")

plot(ts(y[,2]),type="o",ylab="LP")

plot(ts(y[,3]),type="o",ylab="HCT")

par(mfrow=c(1,1))

缺失很多,目的是要补上


5555123

也能有m个时滞(lag),如VAR(m)那样

实际的优点是: 能够适用于各种$A_t$及由矩阵$\Phi$定义的转换,能够拟合更加吝啬的(有较少参数)结构来描述多元时间序列


5555123

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

#-- assumes ss0 has been sourced --#

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

# -- Read Data --

y1=scan("c:/xzwu/2009/berkeley/ts/shumway/mydata/HL.dat")

y2=scan("c:/xzwu/2009/berkeley/ts/shumway/mydata/Folland.dat")

# -- Function to Calculate Likelihood --

Linn=function(para){

cQ=para[1] # sigma_w

cR1=para[2] # 11 element of chol(R)

cR2=para[3] # 22 element of chol(R)

cR12=para[4] # 12 element of chol(R)

cR=matrix(c(cR1,0,cR12,cR2),2) # put the matrix together

kf=Kfilter0(num,y,A,mu0,Sigma0,Phi,cQ,cR)

return(kf$like)

}

# -- Setup --

y=cbind(y1,y2) # y is data matrix (108 x 2)

num=nrow(y)

A=matrix(1,2,1) # A is a 2x1 matrix of 1s - same for all t

mu0=-.35

Sigma0= .01

Phi=1

initpar=c(.1,.1,.1,0) # initial parameter values in order, para[1] to para[4]

# -- Estimation --

#-- view help(optim) for details here...

# the iteration number is written to the screen

# but you have to manually scroll down to see it

#

est=optim(initpar,Linn,NULL,method="BFGS",hessian=TRUE,control=list(trace=1,REPORT=1))

stderr=sqrt(diag(solve(est$hessian)))

#

#-- display summary of estimation

estimate=est$par

u=cbind(estimate,stderr)

rownames(u)=c("sigw","cR11", "cR22", "cR12")

u

#-- Smooth--

#-- first set parameters to their final estimates

cQ=est$par[1]

cR1=est$par[2] # 11 element of chol(R)

cR2=est$par[3] # 22 element of chol(R)

cR12=est$par[4] # 12 element of chol(R)

cR=matrix(c(cR1,0,cR12,cR2),2)

ks=Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,cR)

#-- Plot --

plot.ts(as.vector(ks$xs),lwd=2,ylim=c(-.7,.4), ylab="Temp Deviations")

lines(y1,col="blue",lty="dashed") # color helps here

lines(y2,col="red", lty="dashed")

estimate stderr

sigw 0.04991771 0.01772117

cR11 0.13834951 0.01468864

cR22 0.05793300 0.01009874

cR12 0.04611089 0.01324303


5555123

#Ex6.6

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

#-- assumes ss0 has been sourced --#

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

#-- generate data

set.seed(999)

num=100

N=num+1 # need 101 x's

# generate x(t)=.8x(t-1)+w(t), t=1,...,100

# and x(0) from the stationary distribution:

x = arima.sim(n=N, list(ar = .8, sd=1)) # here you have x0 to x100

v = rnorm(num,0,1) # obs noise

y=ts(x[-1]+v) # observations y[1] ,..., y[100]

#-- initial estimates (method discussed in the text)

u=ts.intersect(y,lag(y,-1),lag(y,-2))

varu=var(u); coru=cor(u)

phi=coru[1,3]/coru[1,2]

q = (1-phi^2)*varu[1,2]/phi

r = varu[1,1] - q/(1-phi^2);

initpar=c(phi,sqrt(q),sqrt(r))

initpar # view the initial estimates

#-- function to evaluate innovations likelihood

Linn=function(para){

phi=para[1]

sigw=para[2] # this is the standard dev

sigv=para[3] # this is the standard dev

mu0=0

Sigma0=(sigw^2)/(1-phi^2)

Sigma0[Sigma0<0] =0 # make sure Sigma0 is never negative

kf = Kfilter0(num,y,1,mu0,Sigma0,phi,sigw,sigv) # run filter under present parameters

return(kf$like) # return -log likelihood

}

#-- do the estimation

#-- view help(optim) for details here...

# the iteration number is written to the screen but you have to manually scroll down to see it

est=optim(initpar, Linn, NULL, method = "BFGS", hessian = TRUE, control=list(trace=1,REPORT=1))

stderr=sqrt(diag(solve(est$hessian))) # standard errors

#-- display summary of estimation

estimate=est$par

u=cbind(estimate,stderr)

rownames(u)=c("phi","sigw","sigv")

u

Run ssall.r

estimate stderr

phi 0.8137623 0.08060636

sigw 0.8507863 0.17528895

sigv 0.8743968 0.14293192


5555123

Global Warming


5555123

#Ex6.10 via BFGS [§6.5]: ex610.txt

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

#-- assumes ss0 has been sourced --#

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!#

# -- Read Data & Make Measurement Matrix

y=ts(scan("c:/xzwu/2009/berkeley/ts/shumway/mydata/jj.dat"),start=1960,freq=4)

num=length(y)

A=cbind(1,1,0,0)

# -- Function to Calculate Likelihood --

Linn=function(para){

phi=para[1]

Phi=diag(0,4); # Phi is 4x4 but only one element is a parameter

Phi[1,1]=phi; Phi[2,]=c(0,-1,-1,-1)

Phi[3,] = c(0,1,0,0); Phi[4,] = c(0,0,1,0)

cQ1=para[2] # sqrt q11

cQ2=para[3] # sqrt q22

cQ=diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2

cR=para[4] # sqrt r11

kf=Kfilter0(num,y,A,mu0,Sigma0,Phi,cQ,cR)

return(kf$like)

}

# -- Initial Parameters --

mu0=c(.7,0,0,0)

Sigma0= diag(.04,4)

initpar=c(1.03,.1,.1,.5) # initial parameters for Phi[1,1], the 2 Q's and R

# -- Estimation --

# the iteration number is printed to the screen but you have to manually scroll to see it

est=optim(initpar,Linn,NULL,method="BFGS",hessian=TRUE,control=list(trace=1,REPORT=1))

stderr=sqrt(diag(solve(est$hessian)))

#-- display summary of estimation

estimate=est$par

u=cbind(estimate,stderr)

rownames(u)=c("Phi11","sigw1","sigw2","sigv")

u

#-- Smooth--

phi=est$par[1]

Phi=diag(0,4);

Phi[1,1]=phi; Phi[2,]=c(0,-1,-1,-1)

Phi[3,] = c(0,1,0,0); Phi[4,] = c(0,0,1,0)

cQ1=est$par[2]

cQ2=est$par[3]

cQ=diag(1,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2 # note lower diag is 2x2 ident for inversions (as a device, they're not used)

cR=est$par[4]

ks=Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,cR)

#-- Plot --

Tsm=ts(ks$xs[1,,],start=1960,freq=4)

Ssm=ts(ks$xs[2,,],start=1960,freq=4)

p1=2*sqrt(ks$Ps[1,1,])

p2=2*sqrt(ks$Ps[2,2,])

par(mfrow=c(3,1))

plot(Tsm, main="Trend Component", ylab="Trend")

lines(Tsm+p1,lty="dashed", col="blue")

lines(Tsm-p1,lty="dashed", col="blue")

plot(Ssm, main="Seasonal Component", ylim=c(-5,4), ylab="Season" )

lines(Ssm+p2,lty="dashed", col="blue")

lines(Ssm-p2,lty="dashed", col="blue")

plot(y, type="p", main="Data (points) and Trend+Season (line)")

lines(Tsm+Ssm)

estimate stderr

Phi11 1.0350847657 0.00253645

sigw1 0.1397255477 0.02155155

sigw2 0.2208782663 0.02376430

sigv 0.0004655672 0.24174702


Dse formulas

Lag operator

dse formulas

ARMA [including ARIMA, VAR: B(L)=I]

input

output

State Space(A linear time-invariant state space representation in innovations form i)

不可观测的状态向量

or(更一般where n_t is the system noise, Q, the system noise matrix, and R the output

(measurement) noise matrix.)


Dse guide pdf dse guide r

注意:看dse-guide.pdf以及 dse-guide.R


5555123

#Shumway Example 6.12

newdat = read.table("c:/xzwu/2009/berkeley/ts/shumway/mydata/Newbold1.txt", comment.char="#")

y = newdat[1:50,2] # quarterly inflation (column 2)

z = newdat[1:50,3] # quarterly interest rates (column 3)

plot(ts(y,start=c(1953,1),frequency=1))

points(ts(y,start=c(1953,1),frequency=1),type='p',pch=15)

lines(ts(z,start=c(1953,1),frequency=1),lty=2)

points(ts(z,start=c(1953,1),frequency=1),type='p',pch=18)

library(dse)

my=TSdata(input= newdat[,3,drop = F],output=newdat[, 2,drop = F])

my=tframed(my, list(start=c(1953,1), frequency=4))

seriesNamesInput(my)="CPI"

seriesNamesOutput(my)="BILLS"

bb=estBlackBox(my,max.lag=3) #lag=?!!!

tfplot(bb)

bb$model #Giving F, G, H, K

bb$estimates #giving $cov, $like, $pred

attributes(bb)#model

#library(dse)

my1=TSdata(output=newdat[,c(3,2),drop = F])

my1=tframed(my1, list(start=c(1953,1), frequency=4))

bb1=estBlackBox(my1,max.lag=3)

tfplot(bb1)

attributes(bb1)#model

bb1$model #Giving F, H, K

bb1$estimates #giving $cov, $like, $pred


5555123

#library(dse)

my2=TSdata(output=newdat[,3,drop = F])

my2=tframed(my2, list(start=c(1953,1), frequency=4))

bb2=estBlackBox(my2,max.lag=3)

tfplot(bb2)

attributes(bb2)#model

bb2$model #Giving F, H, K (one dimention)

bb2$estimates #giving $cov, $like, $pred

bb3=estVARXls(my,max.lag=1)#two variables (one input one output)

bb3$model #Giving A(L), B(L), C(L)

bb3$estimates #giving $cov, $like, $pred

bb4=estVARXls(my1,max.lag=3)#two variables (TWO output)

bb4$model #Giving A(L), B(L), C(L)

bb4$estimates #giving $cov, $like, $pred


Models and computation steps for state space model

Models andComputation Steps for State-Space Model

Sources (for four packages):

Package: FKF, related article (not for R, but having formulas):

http://ideas.repec.org/p/dgr/kubcen/1998141.html (file3826.pdf)

Package: dlm

dlm.pdf (in R project or your R directory)

Package: Shumway:

http://www.stat.pitt.edu/stoffer/tsa2/chap6.htm

Package: dse

In your R directory : dse-guide.pdf and article:

http://www.bank-banque-canada.ca/pgilbert/gil93.pdf


Model shumway notation transition equation and observation equation1

Model (Shumway notation)Transition equation and observation equation

  • xt : p×1 state vector N(m, S)

  • yt : q×1 observation/measurement vector

  • wt : p×1 system noise iid N(0,Q)

  • vt : q×1 observation noise iid N(0,R)

  • ut : r×1 fixed input

  • At : q×p observation/measurement matrix

  • F: p×p ; G: q×r, U: p×r, Q: p×p, R: q×q


Filtering smoothing forecast

Filtering, Smoothing, Forecast

  • s<t: forecasting or prediction

  • s=t: filtering

  • s>t: smoothing


Kalman filter forecast with

Kalman Filter (Forecast ): with

Prediction equations

with

updating equations

where

Kalman gain

The MSE of

is


Kalman smoother with initial

Kalman Smoother: with initial


The lag one covariance smoother

The Lag-one Covariance Smoother

With initial condition

for t=n,n-1,…,2


Three levels of shumway s models

Three “levels” of Shumway’s models


Dse formulas1

dse formulas

ARMA [including ARIMA, VAR: B(L)=I]

State Space

or


Fkf formulas

FKF formulas


Dlm formulas

dlm formulas


Step 1 model specification define model

Step 1: Model Specification(define model)

  • FKF: write function such as arma21ssto create all vectors and matrices etc

  • Shumway: set initial estimates

  • dlm: using functions (and their combination)

    • dlmModARMA , dlmModPoly, dlmModReg , dlmModSeas , dlmModTrig

  • dse: using functions such as ARMA, SS (can transfer to each other) to build an empty model (but)


Step 2 log likelihood function

Step 2: Log Likelihood Function

  • FKF: write function such as objectiveto create log likelihood function

  • Shumway: write function via Kfilter0, Kfilter1, Kfilter2 to create log likelihood function

  • dlm: write function such as buildFun

  • dse: skips this step


Step 3 get estimates from result of step2

Step 3: Get estimates from result of step2

  • FKF: directly using R function optim

  • Shumway directly using R function optim

  • dlm: using function dlmMLE

  • dse: using functions such as estVARXls, estVARXar, estSSfromVARX, estMaxLik, bft, estBlackBox, etc, directly


Step 4 kalman filter filtering smoothing and forecasting

Step 4: Kalman FilterFiltering, Smoothing, and Forecasting

  • FKF: using R function fkf

  • Shumway using functions Kfilter0, Ksmooth0, Kfilter1, Ksmooth1, Kfilter2, Ksmooth2

  • dlm: using function dlmFilter, dlmSmooth, dlmFilterdlmForecast

  • dse: using functions such as forecast, featherForecasts, horizonForecasts, etc

Note: you have to pay attention to the meanings of the output!


Example fkf

Example: FKF

#Simulationdata

ar1 <- 0.6

ar2 <- 0.2

ma1 <- -0.2

sigma <- sqrt(0.2)

## Sample from an ARMA(2, 1) process

a <- arima.sim(model = list(ar = c(ar1, ar2), ma = ma1), n = n, innov = rnorm(n) * sigma)


Example fkf1

Example: FKF

ans <- fkf(a0 = c(0, 0), P0 = diag(c(10, 10)), dt = rbind(0, 0), ct = matrix(0), Tt = matrix(c(ar1, 1, ar2, 0), ncol = 2), Zt = cbind(1, 0), HHt = matrix(c(sigma^2, 0, 0, 0), ncol = 2), GGt = matrix(0), yt = rbind(a))


Example fkf2

Example: FKF


Spss arima

SPSS的实现:ARIMA模型

  • 时间序列的acf和pacf图:可以用选项Graphs-Time Series-Autocorrelations,

  • 然后把变量选入Variables中(对于数据AR1.sav,把时间序列Z选入)。

  • 在Display中(默认地)有选项Autocorrelations和Partial autocorrelations导致acf和pacf图。

  • 人们还经常对残差项绘acf和pacf图。


Spss arima p d q p d q s

SPSS的实现:ARIMA(p,d,q)(P,D,Q)s模型拟合

  • 选择Analyze-Time Series-ARIMA,然后把数据中的时间序列选入Dependent(在数据AR1.sav中,选Z,对数据tssales.sav时选sales,而对数据tsadds2.sav时选sales),对于Independent,仅在使用数据tsadds2.sav时选了adds。

  • 在Model的第一列为ARIMA(p,d,q)(P,D,Q)s模型的前三个参数(p,d,q),第二列(sp,sd,sq)为ARIMA(p,d,q)(P,D,Q)s模型的后三个参数(P,D,Q)。这样只要选定我们所希望尝试的模型参数即可。

  • 周期s由于在定义序列时已经有了(见对话框中注明的Current Periodicity后面的数字),就不用另外输入了。在输出的变量中有误差和拟合(预测)的序列,在输出文件中还有各个参数和一些判别准则等。。


5555123

时间序列的局限性

  • 一个值得研究的时间序列很难能够不受其它时间序列或其它变量的影响。

  • 引入其它变量使得时间序列的数学解决变得困难或者复杂,甚至无法计算。

  • 模型越复杂,假定的条件越多,离实际越远,结果也越不可靠。

  • 有一些其它的模型可以对付这个问题,但解释和内容不同于经典的单变量时间序列方法。

  • 参看dse-guide.pdf以及 dse-guide.R


5555123

公式:指数平滑模型

  • 这些模型中有a,g,d,f为待估计参数,g=0意味着斜率为常数(趋势无变化),而d=0意味着没有季节成分,f和减幅趋势有关;对于时间序列Xt,趋势、光滑后的序列、季节因子和预测的序列分别用Tt、St、It和 表示;另外,p表示周期,et为残差


Linear trend additive seasonality model

指数平滑模型:线性趋势可加季节模型(Linear trend, additive seasonality model)


Linear trend multiplicative seasonality model

指数平滑模型:线性趋势可乘季节模型(Linear trend, multiplicative seasonality model)


Exponential trend additive seasonality model

指数平滑模型:指数趋势可加季节模型(Exponential trend, additive seasonality model)


Exponential trend multiplicative seasonality model

指数平滑模型:指数趋势可乘季节模型(Exponential trend, multiplicative seasonality model)


Damped trend additive seasonality model

指数平滑模型:减幅趋势可加季节模型(Damped trend, additive seasonality model)


Damped trend multiplicative seasonality model

指数平滑模型:减幅趋势可乘季节模型(Damped trend, multiplicative seasonality model)


Arima6

ARIMA模型

  • 平稳时间序列满足的条件:对所有t,E(Zt)=m,而且自协方差函数gts=cov(Xt,Xs)=E(Xt-m)(Xs-m)。 仅仅与差t-s有关,因此可以记gk=gt,t+k=cov(Xt-m)(Xt-m)。对于平稳序列,自相关函数(acf)定义为corr(Zt,Zt+k)= gk/g0。偏相关函数(pacf)定义为corr(Zt,Zt+k|Zt+1,…,Zt+k-1)。

  • 函数acf和pacf的点图可以用来帮助识别平稳过程的ARMA(p,q)模型。AR(p)和MA(q)模型是ARMA(p,q)模型的特例,而ARMA(p,q)模型又是ARIMA(p,d,q)的特例(这里只有趋势,没有季节),而ARIMA(p,d,q)又是既有趋势又有季节成分的ARIMA(p,d,q)(P,D,Q)s模型的特例。


Arima7

ARIMA模型:为了便于描述公式,定义算子


5555123

AR(p)模型

或者,用等价的算子符号,


5555123

MA(q)模型

或者,用等价的算子符号,


Arma p q

ARMA(p,q)模型

或者,用等价的算子符号,


Arima p d q p d q s1

ARIMA(p,d,q)(P,D,Q)s模型

这里 为类似于ARMA(p,q)模型中的算子

只不过是描述季节序列的罢了;它们定义为


5555123

AR(p)模型:

或者,用等价的算子符号,


5555123

AR(p)模型:

或者,用等价的算子符号,


Damped trend multiplicative seasonality model1

指数平滑模型:减幅趋势可乘季节模型(Damped trend, multiplicative seasonality model)


Damped trend multiplicative seasonality model2

指数平滑模型:减幅趋势可乘季节模型(Damped trend, multiplicative seasonality model)


5555123

AR(p)模型:

或者,用等价的算子符号,


  • Login