基于案例学习时间序列 

 时间序列的组成成分 

 系统性部分 

 – 水平 

 – 趋势 

 – 季节性 

 非系统性部分 

 – 噪声/随机扰动 

 时间序列的组成成分 

 加法模型 

 – 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)