R基于案例学习时间序列

基于案例学习时间序列

时间序列的组成成分

系统性部分

– 水平

– 趋势

– 季节性

非系统性部分

– 噪声/随机扰动

时间序列的组成成分

加法模型

– Y = 水平 + 趋势 + 季节性 + 噪声

乘法模型

– Y = 水平 × 趋势 × 季节性 × 噪声

时间序列的可视化

基本方法——时序图  以时间为横坐标，以时间序列相应的取值为纵坐标

局部放大时序图

改变时序图的标度

添加趋势线

剔除季节因素

交互式可视化

library("RODBC")

conn<-odbcConnectExcel2007("Amtrak.xls")

Amtrak<-sqlFetch(conn,"Data")

close(conn)

head(Amtrak)

library(xts)

Amtrak<-xts(Amtrak[,2],order.by=Amtrak[,1]) #创建xts时间序列

plot(Amtrak)

plot(Amtrak['2000/2003'])  #看某区间时间序列图变相放大

Dygraphs JS库——http://dygraphs.com/

– 交互式时间序列图

Dygraphs R包——http://rstudio.github.io/dygraphs/index.html

– 自动绘制时间序列图

– 高度可配置的轴和系列显示

– 丰富的互动功能

– 上下区域显示（如置信带）

library("dygraphs")

lungDeaths<-cbind(mdeaths,fdeaths)  #取数据集

dygraph(lungDeaths)  #绘图 双击返回原图大小,鼠标拖动选取区域放大

#增加时间选择条

dygraph(lungDeaths) %>% dyRangeSelector()

#其他设置

dygraph(lungDeaths) %>% 

  dySeries("mdeaths",label="Male") %>% 

  dySeries("fdeaths",label="Female") %>%

  dyOptions(stackedGraph = TRUE) %>%

  dyRangeSelector(height = 20)

#增加预测局域

hw<-HoltWinters(ldeaths)

predicted<-predict(hw,n.ahead=72,prediction.interval=TRUE)

dygraph(predicted,main="Predicted Lung Deaths (uk)")%>% 

  dyAxis("x",drawGrid=FALSE) %>% 

  dySeries(c("lwr","fit","upr"),label="Deaths") %>% 

  dyOptions(colors=RColorBrewer::brewer.pal(3,"Set1"))  

  

数据预处理

缺失值

– 填补法

– 删除法

相隔时间不相等的时间序列

– 差值法

极端值

用于预测的序列的时间段选择  

###一元线性回归分析-身高与体重

h=c(171,175,159,155,152,158,154,164,168,166,159,164)

w=c(57,64,41,38,35,44,41,51,57,49,47,46)

lm.hw<-lm(h~w)

summary(lm.hw) #摘要

plot(h~w)

abline(lm.hw)

Call:

lm(formula公式 = h ~ w)

Residuals残差:

    Min      1Q  Median      3Q     Max 

-2.9225 -1.3848 -0.1713  1.5498  3.1076 

Coefficients系数:

             Estimate估计值 Std. Error标准差 t value t值 Pr(>|t|)p值显著性检验    

(Intercept截距) 124.36962    3.56286   34.91 8.82e-12 ***

w                 0.79397    0.07391   10.74 8.21e-07 ***

---

Signif有效数字. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error残差标准差: 2.107 on 10 degrees of freedom自由度

Multiple R-squared离差平方和:  0.9203, Adjusted R-squared调整过的离差平方和:  0.9123 

F-statistic F统计量: 115.4 on 1 and 10 DF,  p-value: 8.21e-07

###多元线性回归分析—swiss数据

data(swiss)

head(swiss)

lm.swiss<-lm(Fertility~.,data=swiss) #lm(y~x+k+l) y=ax+bk+cl+d

summary(lm.swiss)

###多元线性回归分析—婚外情数据

library(AER)

data(Affairs)

summary(Affairs)

table(Affairs$affairs) #显示婚外情次数

#数据转化

Affairs$ynaffair[Affairs$affairs>0]<-1  #有婚外情显示1

Affairs$ynaffair[Affairs$affairs==0]<-0 #没有婚外情显示0

Affairs$ynaffair<-factor(Affairs$ynaffair,levels=c(0,1),labels=c("No","Yes")) #给0,1加个标签

table(Affairs$ynaffair) 

#模型构建

fit.full<-glm(ynaffair~.,data=Affairs[,-1],family=binomial())

summary(fit.full)

#模型优化

fit.reduced<- glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())

summary(fit.reduced)

#模型比较

anova(fit.reduced,fit.full,test="Chisq")

#模型参数解释

coef(fit.reduced)

exp(coef(fit.reduced))

#预测

testdata<- data.frame(rating=1:5,age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),

                      religiousness=mean(Affairs$religiousness))

testdata

testdata$prob<-predict(fit.reduced,newdata=testdata,type="response")

testdata  #不同婚姻评价的预测

testdata<-data.frame(rating=mean(Affairs$rating),age=seq(17,57,10),

                     yearsmarried=mean(Affairs$yearsmarried),

                     religiousness=mean(Affairs$religiousness))

testdata   #不同年龄的预测

testdata$prob<-predict(fit.reduced,newdata=testdata,type="response")

testdata

###时间序列的线性趋势

#读取数据

library("RODBC")

conn<-odbcConnectExcel2007("D:/RWork/Amtrak.xls")

Amtrak<-sqlFetch(conn,"Data")

close(conn)

head(Amtrak)

Am<-Amtrak[,c(1,2)]

lm.am<-lm(Ridership~Month,data=Am)

plot(Am)

abline(lm.am)

t<-(1:nrow(Am)) #nrow返回行的个数

lm.am<-lm(Am$Ridership~t)

summary(lm.am)

plot(lm.am)  #残差图均匀分布在0两次最好，最后一个离群值检测

###时间序列的非线性趋势——指数趋势

y<-log(Am$Ridership)

lm.am2<-lm(y~t)

summary(lm.am2)

plot(lm.am2)#par(mfcol=c(2,2))

###时间序列的非线性趋势——多项式趋势

lm.am3=lm(Am$Ridership~t+I(t^2))

summary(lm.am3)

plot(lm.am3)

###随性序列

library("TSA")

data(rwalk)

t=time(rwalk)

model1=lm(rwalk~t)

summary(model1)

plot(rwalk,type='o',ylab='y')

abline(model1) 

###季节性趋势

#季节均值法

# season(tempdub) creates a vector of the month index of the data as a factor 

#season()这个函数创建一个月份向量的数据作为一个因子

data(tempdub)

month.=season(tempdub) # the period sign is included to make the printout from给数据加上时间属性-按月份

# the commands two line below clearer; ditto below.

model2=lm(tempdub~month.-1) # -1 removes the intercept term 消除了截距项

summary(model2)

#有截距的情形，自动将1月的数据丢掉

model3=lm(tempdub~month.) # intercept is automatically included so one month (Jan) is dropped

summary(model3)

#余弦趋势法

har.=harmonic(tempdub,1)

model4=lm(tempdub~har.)

summary(model4)

plot(ts(fitted(model4),freq=12,start=c(1964,1)),ylab='Temperature',type='l',

     ylim=range(c(fitted(model4),tempdub))) 

#参数ylim的设定确保了图象y轴的范围

points(tempdub)