R语言与时间序列学习笔记.docx
《R语言与时间序列学习笔记.docx》由会员分享，可在线阅读，更多相关《R语言与时间序列学习笔记.docx（7页珍藏版）》请在冰豆网上搜索。
R语言与时间序列学习笔记
R语言与时间序列学习笔记
（1）
继续上一次的参数估计话题。
今天分享的是R语言中时间序列的模型初步估计有关内容。
主要有：
时间序列的创建，ARMA模型的建立与模型的参数估计。
一、时间序列的创建
时间序列的创建函数为：
ts（）.函数的参数列表如下：
ts（data=NA,start=1,end=numeric（）,frequency=1,
deltat=1,ts.eps=getOption（"ts.eps"）,class=,names=）
参数说明：
data：
这个必须是一个矩阵，或者向量，再或者数据框frame
Frequency：
这个是时间观测频率数，也就是每个时间单位的数据数目
Start：
时间序列开始值，允许第一个个时间单位出现数据缺失
举例：
ts（matrix（c（NA,NA,NA,1:
31,NA）,byrow=T,5,7）,frequency=7,names=c（"Sun","Mon","Tue","Wen","Thu","Fri","Sat"））
运行上面的代码就可以得到一个日历：
SunMonTueWenThuFriSat
NANANA1234
567891011
12131415161718
19202122232425
262728293031NA
在R语言中本身也有不少数据集，比如统计包中的sunspots，你可以通过函数data（sunspots）来调用它们。
二、一些时间序列模型
这里主要介绍AR，MA，随机游走，余弦曲线趋势，季节趋势等
首先介绍一下AR模型：
AR模型，即自回归（AutoRegressive,AR）模型，数学表达式为：
AR:
y（t）=a1y（t-1）+...any（t-n）+e（t）
其中，e（t）为均值为0，方差为某值的白噪声信号。
那么产生AR模型的数据，我们就有两种方法：
1、调用R中的函数filter（线性滤波器）去产生AR模型；2、根据AR模型的定义自己编写函数
先说第一种方法：
调用R中的函数filter（线性滤波器）去产生AR模型
介绍函数filter的用法如下：
filter（x,filter,method=c（"convolution","recursive"）,
sides=2,circular=FALSE,init）
对于AR
（2）模型x（t）=x（t-1）--0.9x（t-2）+e（t）
w<-rnorm（550）#我们假定白噪声的分布是正态的。
x<-filter（w,filter=c（1,-0.9）,"recursive"）
#方法：
无论是“卷积”或“递归”（可以缩写）。
如果使用移动平均选择“卷积”：
如果“递归”便是选择了自回归。
再说第二种方法：
依据定义自己编程产生AR模型，还是以AR
（2）模型x（t）=x（t-1）--0.9x（t-2）+e（t）为例，可编写函数如下：
w<-rnorm（550）
AR<-function（w）{
x<-w
x[2]=x[1]+w[1]
for（iin3:
550）
x[i]=x[i-1]-0.9*x[i-2]+w[i]
x
}
调用AR（W）即可得到。
如果对相同的随机数，我们可以发现两个产生的时间序列是一致的。
当然对于第二种方法产生的序列需要转换为时间序列格式，用as.ts（）处理。
类似的，我们给出MA，随机游走的模拟：
MA模型：
w<-rnorm（500）
v<-filter（w,sides=2,rep（1,3）/3）
随机游走：
w<-rnorm（200）
x<-cumsum（w）#累计求和，seeexample：
cumsum（1:
!
0）
wd<-w+0.2
xd<-cumsum（wd）
可以做出相应的图形：
再说一下季节性模型：
最简单的季节模型就是一个分段的周期函数。
比如说某地区一年的气温就是一个季节性模型。
利用TSA包里给出的数据tempdub我们可以发现他就是这样的模型
给出验证：
library（TSA）
data（tempdub）
month<-season（tempdub）
model1<-lm（tempdub~month）
summary（model1）
根据R输出的结果：
Call:
lm（formula=tempdub~month）
Residuals:
Min1QMedian3QMax
-8.2750-2.24790.11251.88969.8250
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）16.6080.98716.828<2e-16***
monthFebruary4.0421.3962.8960.00443**
monthMarch15.8671.39611.368<2e-16***
monthApril29.9171.39621.434<2e-16***
monthMay41.4831.39629.721<2e-16***
monthJune50.8921.39636.461<2e-16***
monthJuly55.1081.39639.482<2e-16***
monthAugust52.7251.39637.775<2e-16***
monthSeptember44.4171.39631.822<2e-16***
monthOctober34.3671.39624.622<2e-16***
monthNovember20.0421.39614.359<2e-16***
monthDecember7.0331.3965.0391.51e-06***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
3.419on132degreesoffreedom
MultipleR-squared:
0.9712,AdjustedR-squared:
0.9688
F-statistic:
405.1on11and132DF,p-value:
<2.2e-16
这里2月份系数表明了一月份平均气温与二月份平均气温的差异，以此类推。
在介绍一下一个季节模型：
余弦趋势μ1=βcos（2pi*f*t+φ）
还是考虑上面气温的例子：
验证：
har<-harmonic（tempdub,1）
model2<-lm（tempdub~har）
summary（model2）
看看结果：
Call:
lm（formula=tempdub~har）
Residuals:
Min1QMedian3QMax
-11.1580-2.2756-0.14572.375411.2671
Coefficients:
EstimateStd.ErrortvaluePr（>|t|）
（Intercept）46.26600.3088149.816<2e-16***
harcos（2*pi*t）-26.70790.4367-61.154<2e-16***
harsin（2*pi*t）-2.16970.4367-4.9681.93e-06***
---
Signif.codes:
0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:
3.706on141degreesoffreedom
MultipleR-squared:
0.9639,AdjustedR-squared:
0.9634
F-statistic:
1882on2and141DF,p-value:
<2.2e-16
我们可以作图来看拟合效果：
顺便指出季节模型也可以模拟：
比如μ1=βcos（2pi*f*t+φ）模型可以模拟如下：
t<-1:
500
w<-rnorm（500）
c<-2*cos（2*pi*t/50+0.6*pi+w）
三、自相关与偏自相关
给出自相关的定义：
在信息分析中，通常将自相关函数称之为自协方差方程。
用来描述信息在不同时间
的，信息函数值的相关性。
详情可参见wiki：
http:
//zh.wikipedia.org/wiki/%E8%87%AA%E7%9B%B8%E5%85%B3
我们可以根据定义给出自相关系数（ACF）的算法：
例如数据：
>x<-1:
10
>u<-mean（x）
>v<-var（x）
>sum（（x[1:
9]-u）*（x[2:
10]-u））/（9*v）#延迟1
[1]0.7
>sum（（x[1:
8]-u）*（x[3:
10]-u））/（9*v）#延迟2
[1]0.4121212
>sum（（x[1:
7]-u）*（x[4:
10]-u））/（9*v）#延迟3
[1]0.1484848
在R中也提供了直接计算acf的函数acf（），利用该函数也计算1至3阶的acf，结果如下：
>a<-acf（x,3）
>a
Autocorrelationsofseries‘x’,bylag
0123
1.0000.7000.4120.148
可以看出，是一样的。
利用acf（）可以处理很多阶的acf，以太阳黑子数的数据集做例子：
>data（sunspots）
>acf（sunspots）#给出了相应的图形
>a<-acf（sunspots,6）#为下面做估计做铺垫，列出前6阶的acf
>a
Autocorrelationsofseries‘sunspots’,bylag
0.00000.08330.16670.25000.33330.41670.5000
1.0000.9220.8900.8750.8640.8500.836
偏自相关：
对于一个平稳AR（p）模型，求出滞后k自相关系数p（k）时，实际上得到并不是x（t）与x（t-k）之间单纯的相关关系。
因为x（t）同时还会受到中间k-1个随机变量x（t-1）、x（t-2）、……、x（t-k+1）的影响，而这k-1个随机变量又都和x（t-k）具有相关关系，所以自相关系数p（k）里实际掺杂了其他变量对x（t）与x（t-k）的影响。
为了能单纯测度x（t-k）对x（t）的影响，引进偏自相关系数的概念。
对于平稳时间序列{x（t）}，用数学语言描述就是：
p[（x（t）,x（t-k）]|（x（t-1）,……，x（t-k+1）={E[（x（t）-Ex（t）][x（t-k）-Ex（t-k）]}/E{[x（t-k）-Ex（t-k）]^2}
这就是滞后k偏自相关系数的定义。
总之，偏自相关就是在试图解释在剔除了中间k-1个随机变量x（t-1）、x（t-2）、……、x（t-k+1）的干扰之后，x（t-k）对x（t）影响的相关程度。
在R语言中，使用函数PACF（）可求解
还是使用太阳黑子数的例子：
>b<-pacf（sunspots,6）
>b
Partialautocorrelationsofseries‘sunspots’,bylag
0.08330.16670.25000.33330.41670.5000
0.9220.2720.1890.1350.0640.044
最后，我们利用这两个函数来看看AR（p）,MA（q）的自相关函数与偏自相关函数的截尾性与拖尾性。
利用二中所介绍的方法生成AR
（2），MA
（2）的数据。
AR
（2）模型：
w<-rnorm（550）#我们假定白噪声的分布是正态的。
x<-filter（w,filter=c（1,-0.9）,"recursive"）
MA（3）模型：
w<-rnorm（500）
v<-filter（w,sides=2,rep（1,3）/3）
>qq<-pacf（x,5）
>qq
Partialautocorrelationsofseries‘x’,bylag
12345
0.532-0.861-0.0820.000
可以看出AR
（2）模型的偏自相关函数是截尾的（但由于这个是数据，所以出现pacf只能看出趋势，而不是在2步后直接变为0）
对于MA（3）模型的自相关函数，由于v的第一项与最后一项缺失，不妨截取v的一部分数据，命名为a,有：
>y<-acf（a,5）
>y
Autocorrelationsofseries‘a’,bylag
012345
1.0000.6520.3970.0590.0670.035
也可以看出趋势。
升级会员 