- 外文名
- time series
- 类 型
- 序列
- 作 用
- 描述现象的发展状态和结果
- 要素二
- 指标数值
- 中文名
- 时间序列
- 要素一
- 时间
- 构成要素
- 现象所属的时间等
-
时间序列分析和预测(含实例及代码)
2018-09-17 21:37:34研究时间序列主要目的:进行预测,根据已有的时间序列数据预测未来的变化。 时间序列预测关键:确定已有的时间序列的变化模式,并假定这种模式会延续到未来。 时间序列预测法的基本特点 假设事物发展趋势会...导论
研究时间序列主要目的:进行预测,根据已有的时间序列数据预测未来的变化。
时间序列预测关键:确定已有的时间序列的变化模式,并假定这种模式会延续到未来。
时间序列预测法的基本特点 - 假设事物发展趋势会延伸到未来
- 预测所依据的数据具有不规则性
- 不考虑事物发展之间的因果关系
时间序列数据用于描述现象随时间发展变化的特征。
时间序列分析就其发展历史阶段和所使用的统计分析方法看:传统的时间序列分析和现代时间序列分析。
一、时间序列及其分解
时间序列(time series)是同一现象在不同时间上的相继观察值排列而成的序列。根据观察时间的不同,时间序列中的时间可以是可以是年份、季度、月份或其他任何时间形式。
时间序列:
(1)平稳序列(stationary series)
是基本上不存在趋势的序列,序列中的各观察值基本上在某个固定的水平上波动,在不同时间段波动程度不同,但不存在某种规律,随机波动
(2)非平稳序列(non-stationary series)
是包含趋势、季节性或周期性的序列,只含有其中一种成分,也可能是几种成分的组合。可分为:有趋势序列、有趋势和季节性序列、几种成分混合而成的复合型序列。
趋势(trend):时间序列在长时期内呈现出来的某种持续上升或持续下降的变动,也称长期趋势。时间序列中的趋势可以是线性和非线性。
季节性(seasonality):季节变动(seasonal fluctuation),是时间序列在一年内重复出现的周期波动。销售旺季,销售淡季,旅游旺季、旅游淡季,因季节不同而发生变化。季节,不仅指一年中的四季,其实是指任何一种周期性的变化。含有季节成分的序列可能含有趋势,也可能不含有趋势。
周期性(cyclicity):循环波动(cyclical fluctuation),是时间序列中呈现出来的围绕长期趋势的一种波浪形或振荡式波动。周期性是由商业和经济活动引起的,不同于趋势变动,不是朝着单一方向的持续运动,而是涨落相间的交替波动;不同于季节变动,季节变动有比较固定的规律,且变动周期大多为一年,循环波动则无固定规律,变动周期多在一年以上,且周期长短不一。周期性通常是由经济环境的变化引起。
除此之外,还有偶然性因素对时间序列产生影响,致使时间序列呈现出某种随机波动。时间序列除去趋势、周期性和季节性后的偶然性波动,称为随机性(random),也称不规则波动(irregular variations)。
时间序列的成分可分为4种:趋势(T)、季节性或季节变动(S)、周期性或循环波动(C)、随机性或不规则波动(I)。传统时间序列分析的一项主要内容就是把这些成分从时间序列中分离出来,并将它们之间的关系用一定的数学关系式予以表达,而后分别进行分析。按4种成分对时间序列的影响方式不同,时间序列可分解为多种模型:加法模型(additive model),乘法模型(multiplicative model)。乘法模型:
二、描述性分析
1、图形描述
2、增长率分析
是对现象在不同时间的变化状况所做的描述。由于对比的基期不同,增长率有不同的计算方法。
(1)增长率(growth rate):增长速度,是时间序列中报告期观察值与基期观察值之比减1后的结果,用%表示。由于对比的基期不同,可分为环比增长率和定基增长率。
环比增长率:是报告期观察值与前一时期观察值之比减1,说明现象逐期增长变化的程度;
定基增长率是报告期观察值与某一固定时期观察值之比减1,说明现象在整个观察期内总的增长变化程度。
设增长率为G: 环比增长率 :
定基增长率 :
(2)平均增长率(average rate of increase):平均增长速度,是时间序列中逐期环比值(环比发展速度)的几何平均数减1的结果:
n:环比值的个数
(3)增长率分析中应注意的问题
i: 当时间序列中的观察出现0或负数时,不宜计算增长率。这种序列计算增长率,要么不符合数学公理,要么无法解释其实际意义。可用绝对数进行分析。
ii: 有些情况下,不能单纯就增长率论增长率,注意增长率与绝对水平结合起来。增长率是一个相对值,与对比的基数值的大小有关。这种情况,计算增长1%的绝对值来克服增长率分析的局限性:
增长1%的绝对值表示增长率每增长一个百分点而增加的绝对数量:
增长1%的绝对值=前期水平/100
三、时间序列预测的程序
时间序列分析的主要目的之一是根据已有的历史数据对未来进行预测。时间序列含有不同的成分,如趋势、季节性、周期性和随机性。对于一个具体的时间序列,它可能含有一种成分,也可能同时含有几种成分,含有不同成分的时间序列所用的预测方法是不同的。预测步骤:
第一步:确定时间序列所包含的成分,确定时间序列的类型
第二步:找出适合此类时间序列的预测方法
第三步:对可能的预测方法进行评估,以确定最佳预测方案
第四步:利用最佳预测方案进行预测
1、确定时间序列成分
(1)确定趋势成分
确定趋势成分是否存在,可绘制时间序列的线图,看时间序列是否存在趋势,以及存在趋势是线性还是非线性。
利用回归分析拟合一条趋势线,对回归系数进行显著性检验。回归系数显著,可得出线性趋势显著的结论。
(2)确定季节成分
确定季节成分是否存在,至少需要两年数据,且数据需要按季度、月份、周或天来记录。可绘图,年度折叠时间序列图(folded annual time series plot),需要将每年的数据分开画在图上,横轴只有一年的长度,每年的数据分别对应纵轴。如果时间序列只存在季节成分,年度折叠时间序列图中的折线将会有交叉;如果时间序列既含有季节成分又含有趋势,则年度折叠时间序列图中的折线将不会有交叉,若趋势上升,后面年度的折线将会高于前面年度的折线,若下降,则后面年度的折线将会低于前面年度的折线。
2、选择预测方法
确定时间序列类型后,选择适当的预测方法。利用时间数据进行预测,通常假定过去的变化趋势会延续到未来,这样就可以根据过去已有的形态或模式进行预测。时间序列的预测方法:传统方法:简单平均法、移动平均法、指数平滑法等,现代方法:Box-Jenkins 的自回归模型(ARMA)。
一般来说,任何时间序列都会有不规则成分存在,在商务和管理数据中通常不考虑周期性,只考虑趋势成分和季节成分。
不含趋势和季节成分的时间序列,即平稳时间序列只含随机成分,只要通过平滑可消除随机波动。因此,这类预测方法也称平滑预测方法。
3、预测方法的评估
在选择某种特定的方法进行预测时,需要评价该方法的预测效果或准确性。评价方法是找出预测值与实际值的差距,即预测误差。最优的预测方法就是预测误差达到最小的方法。
预测误差计算方法:平均误差,平均绝对误差、均方误差、平均百分比误差、平均绝对百分比误差。方法的选择取决于预测者的目标、对方法的熟悉程度。
(1)平均误差(mean error):Y:观测值,F:预测值,n预测值个数
由于预测误差的数值可能有正有负,求和的结果就会相互抵消,这种情况下,平均误差可能会低估误差。
(2)平均绝对误差(mean absolute deviation)是将预测误差取绝对值后计算的平均无擦,MAD:
平均绝对误差可避免误差相互抵消的问题,因而可以准确反映实际预测误差的大小。
(3)均方误差(mean square error):通过平方消去误差的正负号后计算的平均误差,MSE:
(4)平均百分比误差和平均绝对百分比误差
ME,MAD,MSE的大小受时间序列数据的水平和计量单位的影响,有时并不能真正反映预测模型的好坏,只有在比较不同模型对同一数据的预测时才有意义。平均百分比误差(mean percentage error,MPE)和平均绝对百分比误差(mean absolute percentage error,MAPE)则不同,它们消除了时间序列数据的水平和计量单位的影响,是反映误差大小的相对值。
4、平稳序列的预测
平稳时间序列只含有随机成分,预测方法:简单平均法、移动平均法、指数平滑法。主要通过对时间序列进行平滑以消除随机波动,又称平滑法。平滑法可用于对时间序列进行短期预测,也可对时间序列进行平滑以描述序列的趋势(线性趋势和非线性趋势)。
(1)简单平均法:根据已有的t期观察值通过简单平均法来预测下一期的数值。设时间序列已有的t期观察值为
,则t+1期的预测值
为:
到了t+1期后,有了t+1期的实际值,t+1期的预测误差
为:
=
t+2期预测值:
简单平均法适合对较为平稳的时间序列进行预测,即当时间序列没有趋势时,用该方法比较好。但如果时间序列有趋势或季节成分,该方法的预测则不够准确。简单平均法将远期的数值和近期的数值看作对未来同等重要。从预测角度,近期的数值比远期的数值对未来有更大的作用,因此简单平均法预测的结果不够准确。
(2)移动平均法(moving average):通过对时间序列逐期递移求得平均数作为预测值的一种预测方法,有简单移动平均法(simple moving average)和加权移动平均法(weighted moving average).
简单移动平均将最近k期数据加以平均,作为下一期的预测值。设移动平均间隔为k(1<k<t),则t期的移动平均值为:
是对时间序列的平滑结果,通过这些平滑值可描述出时间序列的变化形态或趋势。也可以用来预测。
t+1期的简单移动平均预测值为:
t+2期的简单移动平均预测值为:
移动平均法只使用最近k期的数据,在每次计算移动平均值时,移动的间隔都为k,也适合对较为平稳的时间序列进行预测。应用关键是确定合理的移动平均间隔k。对于同一个时间序列,采用不同的移动间隔,预测的准确性是不同的。可通过试验的方法,选择一个使均方误差达到最小的移动间隔。移动间隔小,能快速反映变化,但不能反映变化趋势;移动间隔大,能反映变化趋势,但预测值带有明显的滞后偏差。
移动平均法的基本思想:移动平均可以消除或减少时间序列数据受偶然性因素干扰而产生的随机变动影响,适合短期预测。
(3)指数平滑法(exponential smoothing)是通过对过去的观察值加权平均进行预测,使t+1期的预测值等t期的实际观察值与t期的预测值的加权的平均值。指数平滑法是从移动平均法发展而来,是一种改良的加权平均法,在不舍弃历史数据的前提下,对离预测期较近的历史数据给予较大权数,权数由近到远按指数规律递减,因此称指数平滑。指数平滑有一次指数平滑法、二次指数平滑法、三次指数平滑法等。
一次指数平滑法也称单一指数平滑法(single exponential smoothing),只有一个平滑系数,且观察值离预测时期越久远,权数变得越小。一次指数平滑是将一段时期的预测值与观察值的线性组合作为t+1时期的预测值,预测模型为:
:平滑系数(
)
t+1期的数据是t期的实际观察值与t期的预测值的加权平均。1期的预测值=1期的观察值
2期预测值:
3期预测值:
4期预测值:
对指数平滑法的预测精度,用均方误差来=衡量:
是t期的预测值
加上用
调整的t期预测误差(
)。
使用指数平滑法时, 关键问题是确定一个合适的平滑系数
,不同的
对预测结果产生不同的影响。
=0,预测值仅仅是重复上一期的预测结果;
=1,预测值就是上一期的实际值;
越接近1,模型对时间序列变化的反应就越及时,因为它给当前的实际值赋予了比预测值更大的权数;
越接近0,给当前的预测值赋予了更大的权数,模型对时间序列变化的反应就越慢。
当时间序列有较大随机波动时,选较大
,以便能很快跟上近期的变化;当时间序列比较平稳时,选较小
。
实际应用中,需考虑预测误差,用均方误差衡量预测误差大小。确定
时,可选择几个
进行预测,然后找出预测误差最小的作为最后的
值。
与移动平均法一样,一次指数平滑法可用于对时间序列进行修匀,以消除随机波动,找出序列的变化趋势。
用一次指数平滑法进行预测时,一般
取值不大于0.5,若大于0.5,才能接近实际值,说明序列有某种趋势或波动过大。
阻尼系数
,阻尼系数越小,近期实际值对预测结果的影响越大,反之,越小。阻尼系数是根据时间序列的变化特性来选取。
5、趋势型序列的预测
时间序列的趋势可分为线性趋势和非线性趋势,若这种趋势能够延续到未来,就可利用趋势进行外推预测。有趋势序列的预测方法主要有线性趋势预测、非线性趋势预测和自回归模型预测。
(1) 线性趋势预测
线性趋势(linear trend)是指现象随着时间的推移而呈现稳定增长或下降的线性变化规律。
趋势方程:
:时间序列
的预测值;
是趋势线斜率,表示时间t 变动一个单位,观察值的平均变动数量
(2) 非线性趋势预测
序列中的趋势通常可认为是由于某种固定因素作用同一方向所形成的。若这种因素随时间推移按线性变化,则可对时间序列拟合趋势直线;若呈现出某种非线性趋势(non-linear trend),则需要拟合适当的趋势曲线。
i: 指数曲线(exponential curve):用于描述以几何级数递增或递减的现象,即时间序列的观察值
按指数规律变化,或者时间序列的逐期观察值按一定的增长率增长或衰减。一般的自然增长及大多数经济序列都有指数变化趋势。
趋势方程:
,
为待定系数
若
>1,则增长率随着时间t的增加而增加;若
<1,则增长率随着时间t的增加而降低;若
>0,
<1,则预测值
逐渐降低以到0为极限。
为确定
,
,可采用线性化手段将其化为对数直线形式,两端取对数:
根据最小二乘原理,按直线形式的常数确定方法求得
,
,求出
,
后,再取其反对数,即可得到
,
。
ii: 多阶曲线:
有些现象变化形态复杂,不是按照某种固定的形态变化,而是有升有降,在变化过程中可能有几个拐点。这时就需要拟合多项式函数。当只有一个拐点时,可拟合二项曲线,即抛物线;当有两个拐点时,需要拟合三阶曲线;有k-1个拐点时,需要拟合k阶曲线。
6、复合型序列的分解预测
复合型序列是指含有趋势、季节、周期和随机成分的序列。对这类序列的预测方法是将时间序列的各个因素依次分解出来,然后进行预测。由于周期成分的分析需要有多年的数据,实际中很难得到多年的数据,因此采用的分解模型为:
预测方法有:季节性多元回归模型、季节自回归模型和时间序列分解法预测。
分解法预测步骤:
第一步:确定并分离季节成分。计算季节指数,以确定时间序列中的季节成分。然后将季节成分从时间序列中分离出去,即用每一个时间序列观察值除以相应的季节指数,以消除季节性。
第二步:建立预测模型并进行预测。对消除了季节成分的时间序列建立适当的预测模型,并根据这一模型进行预测。
第三步:计算最后的预测值。用预测值乘以相应的季节指数,得到最终的预测值。
(1)确定并分离季节成分
季节性因素分析是通过季节指数来表示各年的季节成分,以此描述各年的季节变动模式。
i: 计算季节指数(seasonal index)
季节指数刻画了序列在一个年度内各月或各季度的典型季节特征。在乘法模型中,季节指数以其平均数等于100%为条件而构成的,反映了某一月份或季度的数值占全年平均值的大小。若现象的发展没有季节变动,则各期的季节指数应等于100%;若某一月份或季度有明显的季节变化,则各期的季节指数应大于或小于100%。因此,季节变动的程度是根据各季节指数与其平均数(100%)的偏差程度来测定的。
季节指数计算方法较多,移动平均趋势剔除法步骤:
第一步:计算移动平均值(若是季节数据,采用4项移动平均,月份数据则采用12项移动平均),并对其结果进行中心化处理,即将移动平均的结果再进行一次二项移动平均,即得出中心化移动平均值(CMA)。
第二步:计算移动平均的比值,即季节比率,即将序列的各观察值除以相应的中心化移动平均值,然后计算出各比值的季度或月份平均值。
第三步:季节指数调整。由于各季节指数的平均数应应等于1或100%,若根据第二步计算的季节比率的平均值不等于1,则需要进行调整。具体方法:将第二步计算的每个季节比率的平均值除以它们的总平均值。
ii: 分离季节成分
计算出季节指数后,可将各实际观察值分别除以相应的季节指数,将季节成分从时间序列中分离出去:
结果即为季节成分分离后的序列,反映了在没有季节因素影响下时间序列的变化形态。
iii: 建立预测模型并进行预测
7、时序案例分析
https://blog.csdn.net/mengjizhiyou/article/details/104765862
参考:贾俊平《统计学》第六版
《统计学》pdf及课后答案:链接: https://pan.baidu.com/s/1dZPW0smz2cO-67zfn1BFyQ 提取码: ktxq
-
时间序列
2018-09-02 00:33:43时间序列的定义 所谓时间序列就是按照时间的顺序记录的一列有序数据。对时间序列进行观察、研究、找寻他发展变化的规律,预测他将来的走势就是时间序列分析,时间序列分析方法只适用于近期与短期的预测。 相关特征...
时间序列的定义
所谓时间序列就是按照时间的顺序记录的一列有序数据。对时间序列进行观察、研究、找寻他发展变化的规律,预测他将来的走势就是时间序列分析,时间序列分析方法只适用于近期与短期的预测。
相关特征统计量:
-
均值函数序列:反映的是时间序列每时每刻的平均水平
-
方差函数序列:反映的是时间序列围绕其均值做随机波动时平均的波动程度
-
协方差函数和相关系数度量的是两个不同的事件彼此之间的相互影响程序。而自协方差函数和自相关系数(ACF)度量的是同一事件在不同时期之间的相关程度,形象地讲就是度量自己过去的行为对自己现在的影响
-
偏自相关系数:在剔除了中间个随机变量的干扰后,影响的相关变量。
时间序列常用模型:
- 长期趋势:是时间序列在长时期内呈现出来的持续向上或持续向下的变动。
- 季节变动:是时间序列在一年内重复出现的周期性波动。
- 循环波动:是时间序列呈现出得非固定长度的周期性变动。循环波动的周期可能会持续一段时间,但与趋势不同,它不是朝着单一方向的持续变动,而是涨落相同的交替波动。
- 不规则波动:是时间序列中除去趋势、季节变动和周期波动之后的随机波动。不规则波动通常总是夹杂在时间序列中,致使时间序列产生一种波浪形或震荡式的变动。只含有随机波动的序列也称为平稳序列
时间序列分析方法
1、移动平均法
缺点:预测值总是停留在过去的水平上而无法预计会导致将来更高或更低的波动
N越大,修匀的程度也越大,起伏也越小,有利于消除不规则变动的影响,但同时周期变动难于反映出来。一般N值的取值范围:5<=N<=200。当历史序列的基本变化不大且序列中随机变动成分较多时,N的取值应该大一些,否则N的取值应该小一些。选择N值的一个有效方法是,比较若干模型的预测误差,预测误差最小者为好。
一次移动平均法
-
使用最近N期序列值的平均值作为未来各期的预测结果,每期数据求平均时的作用是一样的
-
适用范围:适用于发展趋势变化不大的情况
计算公式为:
预测标准误差为:
加权移动平均法
-
考虑各期数据的重要性,对近期数据给予较大的权重。
-
需要计算总的平均相对误差,修正预测值
-
适用范围:适用于发展趋势变化不大的情况
计算公式为:
二次移动平均法
-
利用移动平均的滞后偏差的规律来建立直线趋势的预测模型,是一种既能反映趋势变化,又可以有效分离出周期变化的方法。
-
适用范围:适用于具有线性趋势的情况
计算公式为:
####2、指数平滑法####指数平滑法实际上是一种特殊的加权移动平均法。预测值是以前观测值的加权和,且对不同的数据给予不同的权,新数据给较大的权,旧数据给较小的权。
缺点:难以确定指数平滑系数,受主观影响较大
一次指数平滑法
- 适用范围:适用于发展趋势变化不大的情况
计算公式为:
即:最新预测值=前一期预测值+前期预测值产生的误差的修正值。平滑系数的确定:
- 平滑系数a起到一个调节器的作用。a值选取得越大,预测值受近期影响越大;a值选取得越小,预测值受远期影响越大。
- 如果时间序列波动不大,比较平稳,则a应取小一点(如0.1~0.5),以减少修正幅度,使预测模型能包含较长时间序列的信息
- 如果时间序列具有迅速且明显的变动倾向,则a应取大一点,如(0.6~0.8)
初始值的确定:
- 当时间序列的数据较多时,取第一期的实际值为初值
- 当时间序列的数据较少时,取最初几期的平均值为初值
二次指数平滑法
- 适用范围:适用于具有线性趋势的情况
计算公式为:
平稳时间序列
一个时间序列,如果均值没有系统的变化(无趋势)、方差没有系统变化,且严格消除了周期性变化,就称之是平稳的。
根据限制条件的严格程度,平稳时间序列可以分为以下两种类型
-
严平稳时间序列:只有当序列的所有统计特征都不会随时间而变化时才能被称为严平稳。是对序列联合分布的要求,以保证序列所有的统计特征都相同
-
宽平稳时间序列:只要求序列二阶平稳,对于高于二阶的维度没有任何要求
平稳时间序列有以下两个统计性质
-
常数均值
-
自协方差函数和自相关函数只依赖于时间的平移长度而与时间的起止点无关
时间序列的预处理
对序列进行纯随机值和平稳性的检验
1、平稳性的检验
对序列的平稳性有两种检验方法,一种是根据时序图和自相关图显示的特征做出判断的图检验方法;一种是构造检验统计量进行假设检验的方法。图检验方法操作简便,但其主观性较强,最好能用假设检验的方法进行辅助判断。
时序图检验
-
根据平稳时间序列方差和均值为常数的特点,平稳序列的时序图应该显示出该序列始终在一个常数值附近随机波动,并且波动的范围由界的特点。
-
如果观察序列的时序图显示出该序列有明显的趋势性或周期性,那它通常不是平稳序列。
如图所示,销量有明显的上升趋势,所以它一定不是平稳序列。
自相关图检验
-
自相关图的一个坐标轴表示延迟时间数,另一个坐标轴表示自相关系数
-
平稳序列通常具有短期自相关性:随着延迟期数的增加,平稳序列的自相关性会很快衰减到0。反之,非平稳序列的自相关系数衰减向0的速度会比较慢。
三角对称性:具有单调趋势的非平稳序列
长期位于零轴的一边:具有单调趋势的非平稳序列
正弦波动规律:具有周期变化规律的非平稳序列
在零轴附近波动:随机性强的平稳时间序列(很快衰减到0)
单位根检验
如果该统计量的P值小于a时,则可以以1-a的置信水平拒绝原假设,认为该序列为平稳序列;否则,接受原假设,认为该序列为非平稳序列。假设a为0.05,此时p值显著大于0.05,则该序列为非平稳序列。
2、纯随机性的检验
纯随机序列又称为白噪声序列。如果序列值之间没有任何的相关性,这种序列我们称为白噪声序列。白噪声序列没有任何分析的价值
- 非平稳序列一定不是白噪声序列
- 平稳序列如果显示出显著的短期相关性,那么该序列一定不是白噪声序列,否则要对其进行纯随机性检验。
如果该统计量的P值小于a时,则可以以1-a的置信水平拒绝原假设,认为该序列为非白噪声序列;否则,接受原假设,认为该序列为纯随机序列。假设a为0.05,此时p值显著小于0.05,则该序列为非白噪声序列。
平稳时间序列分析
某个时间序列经过预处理,被判定为平稳非白噪声序列,就可以利用ARMA模型进行建模
建模步骤:
非平稳时间序列分析
ARIMA模型的实质就是差分运算与ARMA模型的组合。
差分运算:
建模步骤:
-
-
Matlab时间序列分析
2018-11-13 18:53:46在引入时间序列前,先介绍几个matlab函数 matlab中的gallery函数简析 Matlab 中的 gallery 函数是一个测试矩阵生成函数。当我们需要对某些算法进行测试的时候,可以利用gallery函数来生成各种性质的测试矩阵。其用法...文章目录
时间序列分析需要解决的问题
假设我们有这么一段数据(采样自移动公司某段时间的开户客户数目)
下载地址,包含本实验用到的数据/资料以及全部代码,实验报告,引用请告知,关注并留言获取提取码,或者到我的csdn下载进行下载
我们可以发现,在数据中似乎有种明显的趋势规律,或许还存在一定的周期性。但是,我们却不好描述这种趋势或者说是周期,因为数据中存在一定的随机因素模糊了我们的判断。在进行时间序列数据分析时,我们大部分的目的都是想从中找出可以表达的规律以便于我们对未来做预测。
在这里,我们把时间序列数据分解为以下三项:真实数据=趋势项+周期项+随机项 趋势性:例如随时间变化的一次函数/多次函数/幂函数趋势等 周期项:周期规律 随机项:在这里可以理解为噪声,但是这些噪声也是存在规律的。
不难发现,趋势性和周期项都不难描述,难以刻画的是随机项变化。而如果我们真正能预测随机项的变化,那么我们就基本上可以预测整个数据下一步的走势了。
时间序列分析的步骤
- 观测数据(均值,周期等)
- 数据预处理
- 平稳化–去趋势与去周期,剩余随机项
- 用类似于回归/滑动平均的思想来 拟合随机项
1.判断去趋势去周期后的数据是否平稳
2.计算数据的自相关函数和偏相关函数
3.根据自相关/偏相关函数性质决定选用什么模型来拟合随机项
4.模型定阶和拟合参数的求解
5.模型检验- 模型最终实现功能:将周期项/趋势性/随机项结合,是一个关于时间变化的函数,这个函数可以在一定程度上对未来数据进行预测。
如何实现每个步骤
去趋势/去周期
运用差分算子可以去趋势和周期,其中N阶差分用于具有N阶多项式趋势的数据去趋势,N步差分用于具有N步周期的数据去周期。
偏相关/自相关函数的计算
可以运用matlab自带函数,也可以自己编写。
模型定阶
可以参考AIC/BIC等定阶方法
模型检验
通过计算LB统计量和卡方统计量,判断拟合的误差项是否为白噪声序列,如果是,则检验通过。
一个具体的案例 --移动开户数分析
我们有一份移动公司两年间共计700余天的每日开户树,我们想探寻其中的随机因素。
数据观察
变化趋势
图一.原开户数据趋势变化的观察
可以发现,开户数在局部区间的极值变化具有一定的周期性。但是数据太多且有点杂乱。
为了更好的观察数据的规律,拟将数据进行滑动平均处理,在一定程度上更能看出规律性。并且,依据两年的月数天数变化,设置标记每月天数变化的向量,计算每月,每个季度,每年的平均开户数进行观察。计算并绘制得到下图:
图二.开户数据趋势变化周期的观察
上图红色标记为2012年,蓝色标记为2013年。总的看来均有上升的趋势。图二左上角绘制的是2012和2013的开户数滑动平均观察(使用matlab smooth函数,滑动窗口设置为30),容易发现在年这个层面上趋势具有较大的相关性,一年显然可以是一个大周期的刻画。但是这个周期太大了,如果选取这个周期做去季节性,显然会损失大量数据。因此,我们想找到一个更小的周期来刻画,其实我们还可以发现图一中有明显的小起伏,经验证起伏平均间隔为30即一个月,那么一个月是否为一个周期呢?右上图是按月进行平均的统计值绘制,起伏貌似没有什么规律。但是还是可以作为参考。左下角图绘制的是季度的平均开户数变化,可以发现两个季度为一个周期,因为起伏在两个周期中进行。右下角图绘制的是年平均变化图,如我们所料,随着移动互联网的普及,2013年开户平均数比2012年高。
综上观察,我们可以得出以下结论:
1.一个月应该是一个最小周期,两个季度是一个中周期,一年为一个确定的大周期。
2.开户人数有明显的上升趋势。
3.由观测值我们可以认为存在异常值。数据预处理
2.1缺失值处理
数据中不存在缺失值
2.2异常值处理
拟采用三倍标准差即拉以达法则筛选异常值,考虑到时间序列的短期影响性,用其周围值平均代替。去趋势与去周期
进行了简单的观察和数据预处理后,我们开始正式进行时间序列平稳化的操作。考虑到有两年,我们从单独考虑2012,2013,然后集中考虑2012和2013来进行处理和检验。
3.1去趋势
由于数据变化有一定的集中性,且有滑动平均看来近似可以用一次函数拟合。所以采用一次差分的方法去除数据的趋势项。这里采用diff函数实现一次差分。去趋势后,我们进行了如下的观察:图三.周期为一个月的检验图
我们发现相邻极值点间的时间差距十分接近于一个月,这让我们欣喜万分。
3.2去周期
由以上的观察数据和去趋势后数据综合考虑,采用一个月作为周期是个不错的选择。这里自定义d步差分函数对去趋势后的数据进行去周期的处理。
经过3.1和3.2,我们得到如下结果:
图四.2012年去趋势和去周期比较图
图五.2013年去趋势和去周期比较图
图六.2012-2013年去趋势和去周期比较图
肉眼看来,效果还算不错,从均值线看来均值均十分接近于0。为了更客观的看待平稳的效果,采用自相关图还进行平稳性的验证。
自相关性平稳性检验结果如下:
图七.自相关图检验平稳性
从自相关图看来,经过去趋势和去周期后的数据的确是平稳的。 经过以上步骤,我们就成功的将非平稳的时间序列转换为了平稳的时间序列,得到了随机误差项的表示。
截至目前以上代码总结:
%% 计算月平均,季度平均,年平均意图找寻周期规律 % 输入开户数据 % 输出各有用的统计量 function [months_mean,seasons_mean,years_mean,month_starts,month_ends,season_starts,season_ends] = Calu_mean(open_nums) t = 1:length(open_nums); nums = open_nums; month_days = [0,31,29,31,30,31,30,31,31,30,31,30,31,31,28,31,30,31,30,31,31,30,31,30,31]; season_days = [0,sum(month_days(2:4)),sum(month_days(5:7)),sum(month_days(8:10)),sum(month_days(11:13))... ,sum(month_days(14:16)),sum(month_days(17:19)),sum(month_days(20:22)),sum(month_days(23:25))]; % calu the month_mean nums months_mean = zeros(1,24); month_starts = zeros(1,24); month_ends = zeros(1,24); for ii = 1:length(month_days)-1 start_day = sum(month_days(1:ii))+1; end_day = sum(month_days(1:ii+1)); month_starts(ii) = start_day; month_ends(ii) = end_day; months_mean(ii) = mean(nums(start_day:end_day)); end
%% 去趋势和周期函数 % 输入预处理后的时序数据向量vector和观察数据变化人为认定的周期T % 输出去趋势后的向量数据detrend_data和去趋势去周期后的向量数据detrend_deT_data function [detrend_data,detrend_deT_data] = Detrend_plot(vector,T) subplot(3,1,1) plot(1:length(vector),vector,'k') hold on plot(0:length(vector),mean(vector)*ones(1,length(0:length(vector))),'r') text(1,mean(vector),'均值') legend('趋势','均值') title('原开户数据趋势') % 去趋势项 detrend_data = diff(vector,1); subplot(3,1,2) plot(1:length(detrend_data),detrend_data,'k') hold on plot(0:length(detrend_data),mean(detrend_data)*ones(1,length(0:length(detrend_data))),'r') legend('趋势','均值') title('去除趋势项后趋势') % 去季节项 detrend_deT_data = zeros(length(detrend_data),1); for i=1:length(detrend_data) if(i<=T) detrend_deT_data(i) = []; else detrend_deT_data(i) = detrend_data(i)-detrend_data(i-T); end end subplot(3,1,3) plot(1:length(detrend_deT_data),detrend_deT_data,'k') hold on plot(0:length(detrend_deT_data),mean(detrend_deT_data)*ones(1,length(0:length(detrend_deT_data))),'r') legend('趋势','均值') title('去除趋势项后和季节项后趋势') end
%% 主函数脚本,运行该脚本将输出以上所有结果 clearvars %% 读取数据并预处理 open_nums = xlsread('移动通知户开户数.xlsx',1,'B2:B732'); %% 计算月平均,季度平均,年平均并绘图找周期规律,并处理异常值 figure(1),[months_mean,seasons_mean,years_mean,month_starts,month_ends,season_starts,season_ends] = Calu_mean(open_nums); [m,n] = find(abs(open_nums-mean(open_nums))>2*std(open_nums)); open_nums(m) = mean(open_nums); %% 尝试进行去趋势和去周期,考虑到存在多重周期,将年份分开按照月为周期去除周期(其中每月为多少天根据两年分别不同) data_2012 = open_nums(1:366); data_2013 = open_nums(367:end); figure(2),[detrend_data2012,detrend_deT_data2012] = Detrend_plot(data_2012,30);% 2012开户数变化趋势 figure(3),[detrend_data2013,detrend_deT_data2013] = Detrend_plot(data_2013,31);% 2013开户数变化趋势 figure(4),[detrend_data,detrend_deT_data] = Detrend_plot(open_nums,30);% 2012-2013开户数变化趋势 figure(5),subplot(411),autocorr(open_nums),title('原数据自相关图')% 自相关图分析 subplot(412),autocorr(detrend_data),title('2012年随机项自相关图') subplot(413),autocorr(detrend_deT_data2013),title('2013年随机项自相关图') subplot(414),autocorr(detrend_deT_data2013),title('2012-2013年随机项自相关图') clearvars month_starts month_ends season_starts season_ends data_2012 data2013
截至目前,我们编写的代码所完成的工作有:
1.观察数据
2.数据预处理
3.去趋势,去周期
4.计算了偏相关函数和自相关函数我们还需要进行的工作是:
1.判断处理后的随机项是否平稳,判断自相关/偏相关函数的拖尾性/截尾性,以确定采用什么模型来拟合随机项。去均值使该时间序列变为零均值平稳过程。
2.确定模型后,确定采用几阶进行拟合,并求解参数
3.模型检验******* 休息好了吧,我们继续
判断去趋势和去周期后数据的平稳性
图一 处理后的数据 图二 处理后数据的20期自相关系数
自相关系数定义为:
由图二观察知道,当|h|>1时,自相关系数为近似为0.认为自协方差函数(自相关函数)有“截尾性”.可以使用MA(1)模型。为了考虑能否使用ARMA模型,进一步进行讨论。转换为零均值平稳序列
采用mean函数计算该序列的均值,处理如下:
随机变量总体X 的样本的均值估计量:
转换前后序列变化如下:
图三 非零均值平稳序列 图四 零均值平稳序列
3.无偏/有偏自相关函数的估计
自协方差函数(或自相关函数)的样本估计量通常有两种类型:
(3.3.1)
(3.3.2)
关于两种估计量,前人曾经做过比较,G.M.Jenkins 结论:
1.估计量( 3.3.1) 更满意。
2.其中(3.3.1)是自相关函数的有偏估计,(3.3.2)是自相关函数的无偏估计。
3.(3.3.1)式定义的估计量收敛于零的速度更快。
4.由(3.3.1)式定义的样本自协方差函数矩阵(或相关矩阵)是正定的根据两种自相关系数的计算公式,我们可以自己编写函数计算自相关函数如下
其中函数输入为零均值的平稳序列,输出为一组无偏估计的自相关系数,一组为有偏的自相关系数。完成自己的函数编写后,选择计算前面30期自相关系数与matlab自带autocorr进行比较。
自定义计算无偏和有偏两种自相关系数代码:function [unbiased_autocorr,biased_autocorr] = my_autocorr(time_vector) %% 函数定义,计算无偏和有偏两种自相关系数 n = length(time_vector); unbiased_autocorr = zeros(n,1); biased_autocorr = zeros(n,1); % 延迟期数从0到n for h = 0:n omega = 0; % 累加从t=1到n-h for t = 1:n-h omega = omega+time_vector(t)*time_vector(t+h); end unbiased_autocorr(h+1) = omega/(n-h); biased_autocorr(h+1) = omega/n; end % 自相关函数转换为自相关系数 unbiased_autocorr = unbiased_autocorr/unbiased_autocorr(1); biased_autocorr = biased_autocorr/biased_autocorr(1); %% 绘图比较 n_compare = 30; figure() subplot(311) plot(0:n_compare-1,unbiased_autocorr(1:n_compare),'ro','MarkerFaceColor','r','markersize',4) hold on,plot(0:n_compare-1,zeros(1,n_compare),'b') title('无偏自相关系数变化'),xlabel('延迟天数'),ylabel('自相关系数') subplot(312) plot(0:n_compare-1,biased_autocorr(1:n_compare),'ro','MarkerFaceColor','r','markersize',4) hold on,plot(0:n_compare-1,zeros(1,n_compare),'b') title('有偏自相关系数变化'),xlabel('延迟天数'),ylabel('自相关系数') subplot(313) autocorr(time_vector,n_compare),title('内置autocorr结果')
结果如图:
图五 前30期自相关系数比较图
由图,计算结果与内置函数结果十分接近,验证了自己构建的函数的准确性。
图六 无偏/有偏自相关系数的收敛速度比较
可以看出。该图验证了第三个结论,即有偏估计定义的估计量收敛于零的速度更快,且从收敛趋势看来效果更好。
4.偏相关系数的估计和比较
计算偏相关系数有两种方法,直接求解和递推求解。这里考虑递推求解。运用3求解得到的自相关系数一步步递推得到所有的偏相关系数。
4.1 直接求解
(3.4.1)
由最优线性估计的计算公式,应用矩阵按行(列)展开式的性质可得
(3.4.2)
4.2 递推求解
(3.4.3)
递推过程:
图七 递推求解过程
理解以上递推思想后,利用3计算出来的自相关系数,编写函数如下:
其中输入为自相关系数和之前的零均值序列。前者用于递推,后者用于利用matlab自带函数parcorr对自定义的函数计算值进行检验。
代码:function [customed_parcorr] = my_parcorr(time_vector,auto_corr) %% 函数定义 n = length(auto_corr); customed_parcorr = zeros(n,n); customed_parcorr(1,1) = auto_corr(2); for k = 1:n-2 sum1 = 0; sum2 = 0; for j = 1:k sum1 = sum1 + auto_corr(k+2-j)*customed_parcorr(k,j); sum2 = sum2 + auto_corr(j+1)*customed_parcorr(k,j); end customed_parcorr(k+1,k+1) = (auto_corr(k+2)-sum1)/(1-sum2); for j = 1:k customed_parcorr(k+1,j) = customed_parcorr(k,j) - customed_parcorr(k+1,k+1)*customed_parcorr(k,k+1-j); end end %% 绘图比较 a = diag(customed_parcorr); subplot(211) compare_n = 40; parcorr(time_vector,compare_n) subplot(212) plot(1:compare_n,a(1:compare_n),'ro','markerfacecolor','r','markersize',4) hold on,plot(1:compare_n,zeros(1,compare_n),'b') ylim([-0.5,1]) title('自定义偏相关计算值') xlabel('延迟天数'),ylabel('偏相关系数')
结果如图八所示:
图八 自定义偏相关系数计算值与parcorr前30期计算值比较
由比较图看来,计算结果相当。
由图八可知,偏相关系数具有拖尾性。截至目前,除了上面给过的代码外,新的部分计算代码给出如下:
%% 主函数运行脚本 clearvars %% 读取数据并预处理 open_nums = xlsread('移动通知户开户数.xlsx',1,'B2:B732'); %% 计算月平均,季度平均,年平均并绘图找周期规律,并处理异常值 figure(1),[months_mean,seasons_mean,years_mean,month_starts,month_ends,season_starts,season_ends] = Calu_mean(open_nums); [m,n] = find(abs(open_nums-mean(open_nums))>2*std(open_nums)); open_nums(m) = mean(open_nums); %% 尝试进行去趋势和去周期,考虑到存在多重周期,将年份分开按照月为周期去除周期(其中每月为多少天根据两年分别不同) data_2012 = open_nums(1:366); data_2013 = open_nums(367:end); figure(2),[detrend_data2012,detrend_deT_data2012] = Detrend_plot(data_2012,30);% 2012开户数变化趋势 figure(3),[detrend_data2013,detrend_deT_data2013] = Detrend_plot(data_2013,31);% 2013开户数变化趋势 figure(4),[detrend_data,detrend_deT_data] = Detrend_plot(open_nums,30);% 2012-2013开户数变化趋势 figure(5),subplot(411),autocorr(open_nums),title('原数据自相关图')% 自相关图分析 subplot(412),autocorr(detrend_data),title('2012年随机项自相关图') subplot(413),autocorr(detrend_deT_data2013),title('2013年随机项自相关图') subplot(414),autocorr(detrend_deT_data),title('2012-2013年随机项自相关图') clearvars month_starts month_ends season_starts season_ends data_2012 data2013 %% 平稳序列去均值变为零均值平稳序列 zeromean_detrend_deT_data = detrend_deT_data-mean(detrend_deT_data); [unbiased_autocorr,biased_autocorr] = my_autocorr(zeromean_detrend_deT_data); %% 计算偏相关函数 [customed_parcorr] = my_parcorr(zeromean_detrend_deT_data,biased_autocorr);
************** 休息一会,我们继续
处理后数据适用的模型判别
观察平稳序列的自相关函数和偏相关函数,通过是否拖尾/结尾判断使用的模型。
图*1 模型判别流程图
通过上面的分析,自相关系数和自偏相关系数均呈现拖尾性,故使用ARMA模型进行进一步分析。
3.ARMA模型的参数估计和定阶
拟函数法进行求解,求解步骤如下:3.1初定阶数p和q,运用Yule-Walker求解逆函数I 3.2求解MA滑动平均系数 3.3求解AR自回归系数 3.4在求出所有参数后,利用极大似然法残差平方和求方差,算出AIC/BIC指标 3.5重复不走3.1-3.4,在一定大范围内找出极小BIC值对应的阶数p,q
4.参数检验和定阶的进一步修正
通过模型的显著性检验和参数的显著性检验,在参数估计和定阶最小BIC情况下对应的结束p,q附近进行,最终找到最佳的ARMA阶数为ARMA(1,1),参数求解为
经检验残差序列为白噪声,认为模型效果可以接受。
1.适用模型识别
首先,我们需要做的是根据动态数据,从各类模型族中选择与实际过程相吻合的模型;
即处理后的平稳序列应该用什么模型来表示最为合理。经过前期的预处理,去趋势和去周期,得到平稳时间序列的自相关系数,偏相关系数如下图1.1。图1.1 平稳序列的自相关/自协相关系数
不难发现,该平稳序列的自相关系数和自偏相关系数均有拖尾的性质,且自相关拖尾阶数较高分布在32-35阶左右,而自偏相关系数拖尾阶数较低,分布在12-20之间。综上,拟采用ARMA模型对其进行拟合。
2.ARMA模型的参数估计和定阶
确定好适用的模型之后,下一步我们要进行的是对确定的模型类别,定出模型中的阶数p 和 q ,这里采用先大范围搜索,在搜索的结束中求的模型的参数估计——自回归系数和滑动平均系数,并由此求出拟合残差的方差。最后,计算每个阶数下对应的AIC/BIC进行比较,这里适用最小BIC值对应的阶数和参数估计做下一步的分析。2.1 方法的选取
课件4.3关于ARMA模型的参数估计中介绍了两种估计参数的方法,对这两种方法的理解如下:
*ARMA模型参数的矩估计:通过Y-W方程求解出自回归系数,利用MA(q)模型参数的线性迭代法或牛顿—拉裴森算法求得滑动平均系数,最后合并可得拟合模型。这种方法 矩估计方法精度较低,通常只能作为其他迭代参数估计算法的初值。
*模型参数的初估计法—逆函数法:利用可逆域的逆转形式进行求解。这种方法仅涉及到线性方程的求解,计算简单便于实现。
综合上述两种方法,考虑到逆函数法对于低阶或高阶模型都很有效且求解方便。这里采用逆函数法对模型参数进行估计,逆函数法求解步骤如下
逆函数法求解步骤1.令P(*) = max(p,q)+q,用y—w方法估计AR(p*)的 p* 个自回归参数,作为 相应的ARMA( p, q)模型的前 p*个逆函数值. 2.求解线性方程组得到滑动平均系数的估计。 3.将前两步的结果代入课件4.3.6公式,求得自回归系数的估计。
模型参数和求解和定阶
2.2.1模型的定阶指标
*最终预报误差准则(PDE)
其中,方差依赖于k, 其大小反映了模型与数据的拟合程度。第一个因子随 k 增大而增大, 放大了残差方差的不确定性影响.
*最小信息准则(AIC)
AIC准则可应用于ARMA模型及其他统计模型(如多项式回归定阶).能在模型参数极大似然估计的基础上,对于ARMA(p, q)的阶数和相应的参数,同时给出一种最佳估计.
*BIC准则
AIC准则避免了统计检验中由于选取置信度而产生的人为性,为模型定阶带来许多方便. 但AIC方法未给出相容估计.
**注意:AIC阶数估计一般偏高,BIC阶数估计偏低。这里适用BIC准则进行初始阶数的定阶,取BIC值最小下对应的阶数p,q即为初始定阶。
2.2.2 模型的求解
观察之前的自相关系数图和自偏相关系数图,初定ARMA模型初始阶数的选取p和q均在1-20阶内,利用逆函数求解并算的每个p,q下参数对应的BIC值,取最小BIC对应的p,q即为初定的阶数。不同怕p,q下对应的BIC值如下图。
图2.2.2 BIC值的变化
由BIC值的变化我们可以看出,在取ARMA(1,1)时BIC值最小。此时ARMA(1,1)模型建立如下:
定阶的进一步修正-模型的检验
确定好最小BIC准则下的参数估计值和模型阶数p和q后,我们知道——以上各个问题是相互关联的,需整体进行系统化的模型优化. 而根据已有的一组样本数据建立模型,对模型阶数和参数作出判断和估计,是重要的两部分工作.所以最后,我们需要在已经初定的阶数周围小范围内进行搜索和检验,通过模型的有效性检验的参数的显著性检验进一步优化模型,寻找最优的阶数和参数估计值。统计模型只是对生成观测数据真实过程的近似,在模型拟合后,还需要进行模型的有效性检验,及检验拟合模型对序列中信息的提取是否充分。模型的有效性检验即对残差的白噪声检验。
观察残差序列及其自相关系数如下图所示:图3.1 残差序列自相关系数
可以发现序列间没有相关性,不具有拖尾和截尾的性质,可以初步判断改模型是可靠的。为了进一步验证残差序列是否为白噪声序列,下面进行LB统计量的检验结果如下:
延迟期数 LB统计量 卡方统计量
由LB统计量小于对应的卡方统计量,认为该序列可以认为是白噪声序列。模型效果较好,模型检验通过。
到此,我们所有工作已经完毕。如果有什么疑问,欢迎留言,我会一一解答。新部分的代码:
1主函数 main1.m
实现功能:读取处理后的平稳数据进行建模,通过观察自相关/自偏相关系数确定选用ARMA模型,在一定确定的阶数范围内求模型的参数估计,根据最小BIC准则初定p,q的选择,最后通过模型的有效性检验和参数的显著性检验确定最终的参数估计和阶数。% 读取处理后的平稳时间序列 processed_data = csvread('detrend_deT_data.csv'); figure(1) autocorr(processed_data,60) figure(2) parcorr(processed_data,25) data_autocorr = autocorr(processed_data,600); % 自相关系数向量转换为自相关系数矩阵 autocorr_matrix = vec2mat(data_autocorr); max_p = 15; max_q = 20; AICs = ones(max_p,max_q); BICs = ones(max_p,max_q); N = 600; for p=1:20 for q = 1:20 % 计算I P = max(p+q)+q; I = autocorr_matrix(1:P,1:P)\autocorr_matrix(1,2:P+1)'; % 计算MA部分的滑动平均系数 seta = calu_MA_parameter(I,p,q); % 计算AR部分的自回归系数 fine = calu_AR_parameter(seta,I,p,q); % 计算方差 residual_var = calu_var(p,q,fine,seta,N,processed_data); % 计算AIC和BIC并存储 [AICs(p,q),BICs(p,q)] = calu_AIC_BIC(p,q,N,residual_var); end end
% 找出最小BIC下的阶数p,q
[min_BIC_p,min_BIC_q] = find(BICs == min(BICs));
2.MA滑动平均系数的求解函数
[seta]=calu_MA_parameter(I,p,q)
函数输入:逆函数{ Ij }的估计I,给定的自回归和滑动平均阶数p和q
函数输出:滑动平均系数的估计值setafunction [seta]= calu_MA_parameter(I,p,q) P = max(p,q); I_matrix = zeros(q,q); for ii = 1:q I_matrix(ii,:) = I(P+ii-1:-1:P+ii-q); end seta = I_matrix\I(P+1:P+q); end
3.AR自回归系数的求解函数
[fine] = calu_AR_parameter(seta,I,p,q)
函数输入:逆函数{ Ij }的估计I,给定的自回归和滑动平均阶数p和q以及滑动平均系数seta
函数输出:自回归系数的估计值finefunction [fine] = calu_AR_parameter(seta,I,p,q) fine = ones(1,p); for ii = 1:p if ii>q seta(ii) = 0; end fine(ii) = seta(ii)+I(ii); for jj = 1:ii-1 if (ii-jj)>0 fine(ii) = fine(ii)-seta(jj)*I(ii-jj); end end end
4.残差方差的计算函数
[residual_var]= calu_var(p,q,fine,seta,N,processed_data)
函数输入:逆函数{ Ij }的估计I,给定的自回归和滑动平均阶数p和q以及滑动平均系数seta,自回归系数fine,处理后的平稳序列。
函数输出:残差的方差function [residual_var] = calu_var(p,q,fine,seta,N,processed_data) error = zeros(1,N); for ii = max(p,q)+1:N error(ii) = [1,-fine]*processed_data(ii:-1:ii-p)+seta'*error(ii-1:-1:ii-q)'; end residual_var = (error*error')/N; end
5.AIC/BIC指标的计算函数
[AIC,BIC] = calu_AIC_BIC(p,q,N,residual_var)
函数输入:给定的自回归和滑动平均阶数p和q,样本个数N以及残差方差。
函数输出:不同阶数下对应的AIC/BIC指标。function [AIC,BIC] = calu_AIC_BIC(p,q,N,residual_var) AIC = log(residual_var)+2*(p+q+1)/N; BIC = log(residual_var)+log(N)*(p+q+1)/N; End
6.计算LB统计量函数
lb = calu_LB(data, m) 函数输入:残差序列,延迟期数 函数输出:LB卡方统计量 function lb = calu_LB(data, m) N = 25; power_autocorr = 0; data_autocorr = autocorr(data,25); for k=1:m power_autocorr = power_autocorr+data_autocorr(k+1)^2/(N-k); end lb = N*(N+2)*power_autocorr; End
7.自相关系数向量转为矩阵 vec2mat
[autocorr_matrix] = vec2mat(autocorr_vector)
函数输入:自相关系数函数输出:LB卡方统计量function [autocorr_matrix] = vec2mat(autocorr_vector) n = length(autocorr_vector); autocorr_matrix = zeros(n); for ii = 1:n autocorr_matrix(ii,ii:n) = autocorr_vector(1:n-ii+1); end autocorr_matrix = triu(autocorr_matrix)+tril(autocorr_matrix',-1); end
yqol 编写不易,欢迎关注。
-
时间序列预测算法总结
2018-10-18 10:30:48时间序列算法 time series data mining 主要包括decompose(分析数据的各个成分,例如趋势,周期性),prediction(预测未来的值),classification(对有序数据序列的feature提取与分类),clustering(相似数列...时间序列算法
time series data mining 主要包括decompose(分析数据的各个成分,例如趋势,周期性),prediction(预测未来的值),classification(对有序数据序列的feature提取与分类),clustering(相似数列聚类)等。
时间序列的预测
常用的思路:
1、计算平均值
2、exponential smoothing指数衰减
不同的时间点,赋予不同的权重,越接近权重越高
3、snaive:假设已知数据的周期,上一个周期对应的值作为下一个周期的预测值
4、drift:飘移,即用最后一个点的值加上数据的平均趋势
5、Holt-Winters: 三阶指数平滑
Holt-Winters的思想是把数据分解成三个成分:平均水平(level),趋势(trend),周期性(seasonality)。R里面一个简单的函数stl就可以把原始数据进行分解:
一阶Holt—Winters假设数据是stationary的(静态分布),即是普通的指数平滑。
二阶算法假设数据有一个趋势,这个趋势可以是加性的(additive,线性趋势),也可以是乘性的(multiplicative,非线性趋势),只是公式里面一个小小的不同而已。
三阶算法在二阶的假设基础上,多了一个周期性的成分。同样这个周期性成分可以是additive和multiplicative的。 举个例子,如果每个二月的人数都比往年增加1000人,这就是additive;如果每个二月的人数都比往年增加120%,那么就是multiplicative。
性能衡量采用的是RMSE。 当然也可以采用别的metrics:
6、ARIMA: AutoRegressive Integrated Moving Average,ARIMA是两个算法的结合:AR和MA。在ARMA模型中,AR代表自回归,MA代表移动平均。其公式如下:
- 是白噪声,均值为0, C是常数。 ARIMA的前半部分就是Autoregressive:
, 后半部分是moving average:
。 AR实际上就是一个无限脉冲响应滤波器(infinite impulse resopnse), MA是一个有限脉冲响应(finite impulse resopnse),输入是白噪声。
ARIMA里面的I指Integrated(差分)。 ARIMA(p,d,q)就表示p阶AR,d次差分,q阶MA。
为什么要进行差分呢? ARIMA的前提是数据是stationary的,也就是说统计特性(mean,variance,correlation等)不会随着时间窗口的不同而变化。用数学表示就是联合分布相同。
当然很多时候并不符合这个要求,例如这里的airline passenger数据。有很多方式对原始数据进行变换可以使之stationary:
(1) 差分,即Integrated。 例如一阶差分是把原数列每一项减去前一项的值。二阶差分是一阶差分基础上再来一次差分。这是最推荐的做法。
(2)先用某种函数大致拟合原始数据,再用ARIMA处理剩余量。例如,先用一条直线拟合airline passenger的趋势,于是原始数据就变成了每个数据点离这条直线的偏移。再用ARIMA去拟合这些偏移量。
(3)对原始数据取log或者开根号。这对variance不是常数的很有效。
如何做平稳性检验呢?一种方法是肉眼观察,根据肉眼及人的经验判断,不同的人可能会给出不同的结论,还有另一种更客观的统计学方法来辅助判断,也就是“单位根检验”。
对差分后的数据做单位根检验,参考代码如下:
print("单位根检验:\n")
print(ADF(data.diff1.dropna()))
得到:1%、%5、%10不同程度拒绝原假设的统计值和ADF Test result的比较,p-value接近于0,ADF Test result同时小于5%、10%即说明很好地拒绝该假设,也就是差分后的序列是平稳的。
序列平稳后,还需要做白噪声检验
from statsmodels.stats.diagnostic import acorr_ljungbox
acorr_ljungbox(data.diff1.dropna(), lags = [i for i in range(1,12)],boxpierce=True)
结果通过P<α,拒绝原假设,说明差分后的序列是平稳的非白噪声序列,可以进行下一步建模。
模型定阶
经过平稳性检测和白噪声检测,得到平稳的时间序列后,就开始选择合适的ARIMA模型,即ARIMA模型中合适的p,q
第一步我们要先检查平稳时间序列的自相关图和偏自相关图。通过sm.graphics.tsa.plot_acf和sm.graphics.tsa.plot_pacf得到图形
这里就要用到两个很常用的量了: ACF(auto correlation function)和PACF(patial auto correlation function)。对于non-stationary的数据,ACF图不会趋向于0,或者趋向0的速度很慢。
acf(train) ——原始数据
acf(diff(train,lag=1)) ——一阶差分
acf(diff(diff(train,lag=7))) ——去除周期性的一阶差分
确保stationary之后,下面就要确定p和q的值了。定这两个值还是要看ACF和PACF:
AR(p)模型,PACF会在lag=p时截尾,也就是,PACF图中的值落入宽带区域中。
MA(q)模型,ACF会在lag=q时截尾,同理,ACF图中的值落入宽带区域中。
确定好p和q之后,就可以调用R里面的arime函数了。
ARIMA更多表示为 ARIMA(p,d,q)(P,D,Q)[m] 的形式,其中m指周期(例如7表示按周),p,d,q就是前面提的内容,P,D,Q是在周期性方面对应的p,d,q含义。
R里面有两个很强大的函数: ets 和 auto.arima。 用户什么都不需要做,这两个函数会自动挑选一个最恰当的算法去分析数据。
在R中各个算法的效果如下:
模型优化
其中L是在该模型下的最大似然,n是数据数量,k是模型的变量个数。
python代码如下:
arma_mod20 = sm.tsa.ARIMA(data["xt"],(1,1,0)).fit()
arma_mod30 = sm.tsa.ARIMA(data["xt"],(0,1,1)).fit()
arma_mod40 = sm.tsa.ARIMA(data["xt"],(1,1,0)).fit()
values = [[arma_mod20.aic,arma_mod20.bic,arma_mod20.hqic],[arma_mod30.aic,arma_mod30.bic,arma_mod30.hqic],[arma_mod40.aic,arma_mod40.bic,arma_mod40.hqic]]
df = pd.DataFrame(values,index=["AR(1,1,0)","MA(0,1,1)","ARMA(1,1,0)"],columns=["AIC","BIC","hqic"])
df
构造这些统计量所遵循的统计思想是一致的,就是在考虑拟合残差的同时,依自变量个数施加“惩罚”。但要注意的是,这些准则不能说明某一个模型的精确度,也即是说,对于三个模型A,B,C,我们能够判断出B模型是最好的,但不能保证B模型能够很好地刻画数据。可以用AIC(mod1)来检验与比较,值越小越好。
参数估计
from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(data["xt"], order=(0,1,1))
result = model.fit()
print(result.summary())
模型检验
P<α,拒绝原假设,认为该参数显著非零MA(2)模型拟合该序列,残差序列已实现白噪声
resid = result.resid#残差
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)
fig = qqplot(resid, line='q', ax=ax, fit=True)
qq图显示,我们看到红色的KDE线与N(0,1)平行,这是残留物正太分布的良好指标,说明残差序列是白噪声序列,模型的信息的提取充分,当让大家也可以使用前面介绍的检验白噪声的方法LB统计量来检验
模型预测
pred = result.predict('1988', '1990',dynamic=True, typ='levels')
print (pred)
plt.figure(figsize=(12, 8))
plt.xticks(rotation=45)
plt.plot(pred)
plt.plot(data.xt)
plt.show()
预测效果评价:accuracy()函数包含了很多评价指标,也可以自定义。
7、ARFIMA(p,d,q):分数积分自回归移动平均模型,时间序列长期记忆
关键步骤在于对序列进行分数阶差分,得到满足零均值ARMA过程的序列,再进行参数估计。
案例:
我们有一个时间序列的数据 "D1_data.txt"
1.拿到这个时间序列你先目测,上图通过目测可以认为是一个平稳的时间序列。
R步骤
a <- read.csv(“D:/D1_data.txt”) #括号中的东西取决于你的数据在电脑上存放的路径
aa <-ts(a) #命名该事件序列为aa
plot(aa) #画出时间序列aa如上图
2、通过ACF和PACF也就是自相关系数和偏自相关系数来判断是否稳定
由于ACF在lag=1之后便落入置信区间,PACF在lag=3之后落入置信区间,我们认为该时间序列稳定。
R步骤
layout(1:2) #一次显示两张图
acf(aa) #显示自相关系数
pacf(aa) #显示偏自相关系数
3、模型拟合
选择ma(1)或者arma(3,1)
用arima命令来拟合,并用AIC来看哪个模型更适合,AIC数值越小越好。
R步骤
ma1 <- arima(aa,order = c(0,0,1))
arma13 <- arima(aa,order = c(3,0,1))
AIC(ma1)
AIC(arma13)
#AIC结果:
AIC(ma1) [1] 939.6636
AIC(arma13) [1] 943.8848
可以看出ma1优于arma13.
4、看拟合好的模型的残差是否为白噪声
R步骤
acf(resid(ma1)) #ma1的残差的自相关系数
注意事项:
对于这几个模型的选择绝对不是IC谁高就选谁。
正确的流程:1)根据acf和pacf的特点进行选择。比如题主的图acf缓慢下行,而pacf只有一个lag significant,这样对应的model是AR,相反如果pacf decays而acf cuts off,你就该去选MA。每个model都有对应的acf,pacf的特点,根据对应的特点细节去选model,选lag,增加sessional ar ma,验证残差等等。
2)在1)的基础上选出的可能性的model中,用IC和残差的表现来选一个最好的,这时候才用到IC。不谈1)来谈IC没有意义。其中还有很多需要注意的点建议题主自己搜索学习一下。
有trend代表一定不stationary,但这并不影响模型,在去除趋势项后,一个时间序列依然可以稳态,这个叫做trend stationary,因此即使有trend,依然可以使用ar等模型来model。
存在趋势的序列都是非平稳的,AR等一系列模型是必须建立在平稳的基础上才有意义…一般时间序列建模的流程是:去除确定性因素(趋势还有季节性),然后对剩下的随机因素进行平稳性检验,检验通过之后进行arima建模,具体的阶数你可以用acf,pacf来确定,比较方便的是R里面的auto.arima直接定阶…然后进行模型的拟合与预测…最后就是对残差进行arch效应检验,如果有arch效应,就进一步用波动率模型进行建模…总的来说就是确定性因素的分解—随机因素的均值建模—波动率建模。
疑问:
如何判断一个序列是否平稳?
1、均值 ,是与时间t 无关的常数。
2、方差 ,是与时间t 无关的常数。这个特性叫做方差齐性。
3、协方差 ,只与时期间隔k有关,与时间t 无关的常数。
如果一个序列不平稳,就要用方法使它变平稳:如消除长期趋势、差分化
检查平稳性的公式,
1、引入Rho系数:X(t) = Rho * X(t-1) + Er(t)
2、Dickey-Fuller检验是测试一个自回归模型是否存在单位根。X(t) - X(t-1) = (Rho - 1) X(t - 1) + Er(t),要测试如果Rho–1=0是否差异显著。如果零假设不成立,我们将得到一个平稳时间序列。
如何判断一个序列符合AR还是MA?
自回归AR案例:
例如,x(t)代表一个城市在某一天的果汁的销售量。在冬天,极少的供应商进果汁。突然有一天,温度上升了,果汁的需求猛增到1000.然而过了几天,气温有下降了。但是众所周知,人们在热天会喝果汁,这些人会有50%在冷天仍然喝果汁。在接下来的几天,这个比例降到了25%(50%的50%),然后几天后逐渐降到一个很小的数。
移动平均MA案例:
一个公司生成某种类型的包,这个很容易理解。作为一个竞争的市场,包的销售量从零开始增加的。所以,有一天他做了一个实验,设计并制作了不同的包,这种包并不会被随时购买。因此,假设市场上总需求是1000个这种包。在某一天,这个包的需求特别高,很快库存快要完了。这天结束了还有100个包没卖掉。我们把这个误差成为时间点误差。接下来的几天仍有几个客户购买这种包。
AR模型与MA模型的不同
AR与MA模型的主要不同在于时间序列对象在不同时间点的相关性。
MA模型用过去各个时期的随机干扰或预测误差的线性组合来表达当前预测值。当n>某一个值时,x(t)与x(t-n)的相关性总为0.
AM模型仅通过时间序列变量的自身历史观测值来反映有关因素对预测目标的影响和作用,步骤模型变量相对独立的假设条件约束,所构成的模型可以消除普通回退预测方法中由于自变量选择、多重共线性等造成的困难。即AM模型中x(t)与x(t-1)的相关性随着时间的推移变得越来越小。
用ACF和PACF来区分AR与MA
时间序列x(t)滞后k阶的样本自相关系数(ACF)和滞后k期的情况下样本偏自相关系数(PACF)。
AR模型的ACF和PACF:
通过计算证明可知:
- AR的ACF为拖尾序列,即无论滞后期k取多大,ACF的计算值均与其1到p阶滞后的自相关函数有关。
- AR的PACF为截尾序列,即当滞后期k>p时PACF=0的现象。
上图蓝线显示值与0具有显著的差异。很显然上面PACF图显示截尾于第二个滞后,这意味这是一个AR(2)过程。
MA模型的ACF和PACF:
- MA的ACF为截尾序列,即当滞后期k>p时PACF=0的现象。
- AR的PACF为拖尾序列,即无论滞后期k取多大,ACF的计算值均与其1到p阶滞后的自相关函数有关。
很显然,上面ACF图截尾于第二个滞后,这以为这是一个MA(2)过程
如何让时间序列变得平稳?
1 消除趋势:这里我们简单的删除时间序列中的趋势成分
2 差分:这个技术常常用来消除非平稳性。这里我们是对序列的差分的结果建立模型而不是真正的序列。
x(t) – x(t-1) = ARMA (p , q)
差分后,ARIMA对应为:p:AR d:I q:MA
3 季节性:季节性直接被纳入ARIMA模型中。
如何确定pdq的值?
1、参数p,q可以使用ACF和PACF图发现。
2、如果相关系数ACF和偏相关系数PACF逐渐减小,这表明我们需要进行时间序列平稳并引入d参数。
选择模型时,选择AIC和BIC最小的(p,d,q)组合。
进行时间序列预测的步骤
- 是白噪声,均值为0, C是常数。 ARIMA的前半部分就是Autoregressive:
-
时间序列的7种预测模型
2018-12-23 20:02:00时间序列问题比较常见,比如股市,工业生产指标等。 1 朴素估计 使用最后一个时间点的值估测后面一段时间段的值。 2 简单平均 4 滑动窗平均 使用之前一定大小时间段的平均值作为这个时间点的值。 ... -
时间序列分析这件小事(六)--非平稳时间序列与差分
2016-12-04 21:31:29之前我们说明了怎么样的时间序列是序列平稳的,但是世界并不是那么美好,很多时间序列都不是平稳序列,所以这里就要求我们做一些处理了。 首先我们来看一下非平稳时间序列长什么样。在AR模型中,只要自回归系数都... -
时间序列matlab的实现
2019-08-23 11:14:45时间序列预测法其实是一种回归预测方法,属于定量预测,其基本原理是:一方面承认事物发展的延续性,运用过去的时间序列数据进行统计分析,推测出事物的发展趋势;另一方面充分考虑到由于偶然因素影响而产生的随机性... -
时间序列模型 (七): 时间序列建模的基本步骤
2019-04-22 12:21:50时间序列建模的基本步骤 时间序列建模的基本步骤 习题 时间序列经典教材推荐 时间序列模型 (一):模型概述 时间序列模型 (二):移动平均法 时间序列模型 (三):指数平滑法 时间序列模型 (四):差分... -
简单粗暴LSTM:LSTM进行时间序列预测
2018-12-14 21:12:48使用LSTM进行时间序列预测 欢迎使用Markdown编辑器 你好! 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章,了解一下Markdown的基本语法知识。 新... -
时间序列分析
2017-03-22 17:04:51http://blog.csdn.net/pipisorry/article/details/62053938时间序列简介时间序列是时间间隔不变的情况下收集的时间点集合。这些集合被分析用来了解长期发展趋势,为了预测未来或者表现分析的其他形式。但是什么时间... -
时间序列分析之ADF检验
2019-02-06 18:59:09在使用很多时间序列模型的时候,如 ARMA、ARIMA,都会要求时间序列是平稳的,所以一般在研究一段时间序列的时候,第一步都需要进行平稳性检验,除了用肉眼检测的方法,另外比较常用的严格的统计检验方法就是ADF检验... -
算法模型---时间序列模型
2018-01-16 09:04:581、时间序列时间序列是时间间隔不变的情况下收集的不同时间点数据集合,这些集合被分析用来了解长期发展趋势及为了预测未来。 时间序列与常见的回归问题的不同点在于: 1、时间序列是跟时间有关的;而线性回归模型... -
【时间序列】时间序列分解总结
2020-10-22 18:22:00作者| 追光者研究| 机器学习与时间序列出品 | AI蜗牛车1.朴素分解一个时间通常由长期趋势,季节变动,循环波动,不规则波动几部分组成长期趋势指现象在较长时期内持续发展变化的一种趋... -
pandas时间序列
2018-07-25 22:12:35时间序列(time series)数据是一种重要的结构化数据形式,。在多个时间点观察或测量到的任何时间都可以形成一段时间序列。很多时间, 时间序列是固定频率的, 也就是说, 数据点是根据某种规律定期出现的(比如每15... -
时间序列相似性度量-DTW
2019-01-25 16:44:59最近项目中遇到求解时间序列相似性问题,这里序列也可以看成向量。在传统算法中,可以用余弦相似度和pearson相关系数来描述两个序列的相似度。但是时间序列比较特殊,可能存在两个问题: 两段时间序列长度不同。... -
Python时间序列分析视频教程
2017-07-21 13:53:07Python时间序列分析视频培训教程,该课程以Python为核心工具,使用Pandas库进行时间序列的预处理与分析。详解时间序列中常用模型ARIMA原理以及其参数选择。对时间序列平稳性以及模型评估方法展开分析讨论,基于真实... -
用python做时间序列预测三:时间序列分解
2020-06-02 16:37:00分解的使用场景有很多,比如当我们需要计算该时间序列是否具有季节性,或者我们要去除该时间序列的趋势和季节性,让时间序列变得平稳时都会用到时间序列分解。 加法和乘法时间序列 时间序列的各个观测值可以是以上... -
【时间序列】时间序列分析基本方法和实例
2020-04-28 18:00:36目录1 数据相关2 时间序列中的模型(Patterns)3 如何分解时间序列中的各个成分4 平稳与不平稳时间序列4.1 这些数据有什么明显的特点?4.2 为什么要在预测前把序列变成平稳的?4.3 如何对平稳性进行测试4.4 白噪音和... -
时间序列分类算法之时间序列森林(TSF)
2020-04-15 19:57:02时间序列森林(Time Series Forest, TSF)模型将时间序列转化为子序列的均值、方差和斜率等统计特征,并使用随机森林进行分类。TSF通过使用随机森林方法(以每个间隔的统计信息作为特征)来克服间隔特征空间巨大的... -
时间序列与非时间序列的异常检测
2018-04-28 10:51:23非时间序列:序列对应的事件发生的时间无明确先后顺序关系时间序列:序列事件的发生要考虑先后顺序关系相对来说,非时间序列的异常检测就是散点的检测,只需要使用距离度量进行聚类、分类等操作,发现事件中的离群点... -
时间序列实例
2019-02-19 14:28:54时间序列实例 junjun 2016年2月12日 Rmarkdown脚本及数据集:http://pan.baidu.com/s/1gekA3AV 实例一、使用ARIMA模型对裙子长度预测 ARIMA 模型为平稳时间序列定义的。 因此, 如果你从一个非... -
时间序列预测
2017-09-03 23:59:57时间序列预测相关概念时间序列 时间序列中连续的、不同时刻的随机变量,他们彼此之间都有一定的相关性 按照时间的顺序把事件变化发展的过程记录下来就构成了一个时间序列 时间序列预测 对时间序列进行观察、研究... -
时间序列分类
2018-08-19 10:15:16时间序列分类比较麻烦是因为我们用于模型训练的数据的每条样本一般是一个特征向量x对应一个y的形式,而时间序列的大量的信息藏在它的结构中,不仅仅体现在数值上。没意识到这一点的话,我们提取的特征可能就没有什么... -
时间序列分析之协整检验
2019-02-07 14:18:02平稳性是进行时间序列分析的一个很重要的前提,很多模型都是基于平稳下进行的,而现实中,很多时间序列都是非平稳的,所以协整是从分析时间序列的非平稳性入手的。 协整的内容是: 设序列是 d 阶单整的,记为,... -
横截面数据、时间序列数据、面板数据
2018-03-20 15:12:40面板数据(Panel Data)是将“截面数据”和“时间序列数据”综合起来的一种数据类型。具有“横截面”和“时间序列”两个维度,当这类数据按两个维度进行排列时,数据都排在一个平面上,与排在一条线上的一维数据有着... -
时间序列预测,非季节性ARIMA及季节性SARIMA
2019-03-24 21:55:00Python 3中使用ARIMA进行时间序列预测的指南 在本教程中,我们将提供可靠的时间序列预测。我们将首先介绍和讨论自相关,平稳性和季节性的概念,并继续应用最常用的时间序列预测方法之一,称为ARIMA。 介绍 时间...
-
C#学生管理系统源码
-
电商设计专业思维
-
【python基础】递归函数
-
【数据分析-随到随学】数据分析建模和预测
-
MFC开发简单聊天程序
-
前端架构师-速成
-
三维地图GIS大数据可视化
-
冷作工工艺与技能训练(第二版)习题册参考资料(答案)-A02-1134.pdf
-
工程力学(第六版)习题册参考答案-A02-3600.pdf
-
焊工工工艺学(第四版)习题册参考答案-A02-0786.pdf
-
ajax
-
Infineon-HV_Floating_MOS_Gate_Drivers-ApplicationNotes-v01_00-CN.pdf
-
焊工工工艺学(第四版)习题册参考答案-A02-0786.pdf
-
TensorFlow cnn 图片分类源码 用这个 足以应付一切
-
Mybatis-plus argument type mismatch
-
Tomcat8+JDK8安装与配置
-
旅游英语(第二版)参考资料(答案)-+A32-2390.pdf
-
店音,连锁店背景音乐管理软件播放器
-
解决pycharm安装torch库失败的问题
-
VScode快速生成模板