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)