時間序列分析R code 參考

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

pm2.5 時間序列分析 R code 參考 (成果在附件中)

Pm2.5模型建立

#-------------------------------------------------------------------

##建立時間序列資料

library('TSA')

pm2.5=read.csv("C:/pm2.52007-2016(mug.m3).csv",header=F)

pm2.5ts=ts(pm2.5[,1],start=c(2007,1),end=c(2015,12),frequency=12)

win.graph(width=4.875, height=3,pointsize=8)

plot(pm2.5ts,ylab='pm2.5',main='pm2.52007-2015 monthly average')

month=c("J","F","M","A","M","J","J","A","S","O","N","D")

points(window(pm2.5ts),pch=month)

#-------------------------------------------------------------------

##0.穩定模型

BoxCox.ar(pm2.5ts, lambda = seq(0, 1,0.02))

pm2.5ts1=sqrt(pm2.5ts)

win.graph(width=4.875,height=3,pointsize=8)

plot(pm2.5ts1,ylab='sqrt(pm2.5)',main='sqrt(pm2.5)2007-2015 monthly average',ylim=c(0,10))

month=c("J","F","M","A","M","J","J","A","S","O","N","D")

points(window(pm2.5ts1),pch=month)

acf(pm2.5ts1,lag.max=36,

main=expression(Sample~~ACF~~of~~pm2.5ts1~~Levels))

trend2 = lm(pm2.5ts1~time(pm2.5ts1))

summary(trend2)

yt =ts(rstudent(trend2),start=c(2007,1),frequency=12)

plot(yt,type='o',main='de-trended timeseries')

month=c("J","F","M","A","M","J","J","A","S","O","N","D")

points(window(yt),pch=month)

#-------------------------------------------------------------------

##a.判定模型 linear trend +SARIMA(0,0,2)*(0,1,1)

par(mfrow=c(1,1))

acf(yt,lag.max=36,

main='de=trended sqrt(pm2.5) 2007-2015monthly average')

pacf(yt)

par(mfrow=c(1,2))

acf(diff(yt,lag=12),lag.max=36,ci.type='ma',

main='sqrt(pm2.5) 2007-2015 monthlyaverage')

pacf(diff(yt,lag=12),lag.max=36,

main='sqrt(pm2.5) 2007-2015 monthly average')

eacf(diff(yt,lag=12))

#-------------------------------------------------------------------

##b.參數估計

m2.pm2.5=arima(yt,order=c(0,0,2),seasonal=list(order=c(0,1,1),period=12))

m2.pm2.5

#-------------------------------------------------------------------

##c.殘差分析

res2 = residuals(m2.pm2.5)

plot(res2)

par(mfrow=c(3,1))

plot(res2)

hist(res2)

qqnorm(res2); qqline(res2)

t.test(res2)

shapiro.test(res2)

par(mfrow=c(1,1))

acf(res2,lag = 36)

Box.test(res2, lag = 23,type ="Ljung")

#-------------------------------------------------------------------

##d.其他候選模型

##SARIMA(0,1,1)*(0,1,1)

pm2.5ts2=diff(pm2.5ts1)

acf(pm2.5ts2,lag.max=36,

main=expression(Sample~~ACF~~of~~pm2.5ts2~~Levels))

pm2.5ts3=diff(pm2.5ts2,lag=12)

acf(pm2.5ts3,lag.max=36,

main=expression(Sample~~ACF~~of~~pm2.5ts2~~Levels))

pacf(pm2.5ts3,lag.max=36,

main=expression(Sample~~ACF~~of~~pm2.5ts2~~Levels))

eacf(pm2.5ts3)

m1.pm2.5=arima(pm2.5ts1,order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12))

m1.pm2.5

m1.pm2.5=arima(pm2.5ts1,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))

m1.pm2.5

res1 = residuals(m1.pm2.5)

par(mfrow=c(3,1))

plot(res1)

hist(res1)

qqnorm(res1); qqline(res1)

t.test(res1)

shapiro.test(res1)

par(mfrow=c(1,1))

acf(res1,lag = 36)

Box.test(res1, lag = 15,type ="Ljung")

Box.test(res1, lag = 23,type ="Ljung")

#-------------------------------------------------------------------

##預測

par(mfrow=c(1,1))

pm2.511=pm2.5[,1]

testpm2.5 =ts(pm2.511[(9*12+1):(10*12)],start=c(2016,1),frequency=12)

testpm2.5=(testpm2.5)^(1/2)

m2.pm2.51=arima(pm2.5ts1,order=c(1,0,0),seasonal=list(order=c(0,1,1),period=12)

,xreg=as.matrix(model.matrix(~time(pm2.5ts1)))[,-1])

newtrend =time(pm2.5ts1)[length(pm2.5ts1)]+(1:12)*deltat(pm2.5ts1)

plot(m2.pm2.51,n1=c(2007,1),n.ahead=12,newxreg= as.matrix(model.matrix(~newtrend))[,-1]

,col='red',xlab='Year',type='o',ylab='sqrt(pm2.5)',

main=expression(Long~~Term~~Forecasts~~'for'~~the~~pm2.5~~Model))

points(testpm2.5,col = "blue",pch=3)

降雨量模型建立

#-------------------------------------------------------------------

##建立時間序列資料

rain=read.csv("C:/rain2007-2016(mm3).csv",header=F)

raints=ts(rain[,1],start=c(2007,1),end=c(2015,12),frequency=12)

win.graph(width=4.875,height=3,pointsize=8)

plot(raints,ylab='CO2',main='rain 2007-2015monthly average')

month=c("J","F","M","A","M","J","J","A","S","O","N","D")

points(window(raints),pch=month)

#-------------------------------------------------------------------

##穩定模型

BoxCox.ar(raints, lambda = seq(0, 1, 0.02))

raints1=(raints)^(1/4)

plot(raints1,ylab='CO2',main='rain2007-2015 monthly average')

month=c("J","F","M","A","M","J","J","A","S","O","N","D")

points(window(raints1),pch=month)

month.=season(raints1)

trend1 = lm(raints~time(raints1)+month.-1)

summary(trend1)

model = lm(raints1~month.-1)

summary(model)

plot(y=rstudent(model),x=as.vector(time(raints1)),ylab='De-trendseries',

xlab='Time',type='o')

yt = rstudent(model)

#-------------------------------------------------------------------

##a.模型判定

par(mfrow=c(2,1))

acf(yt,lag = 36)

pacf(yt,lag = 36)

eacf(yt)

#-------------------------------------------------------------------

##b.參數估計

model = lm(raints1~month.)

summary(model)

model = lm(raints1~month.-1)

summary(model)

#-------------------------------------------------------------------

##c.殘差分析

res1 = residuals(model)

par(mfrow=c(3,1))

plot(res1)

hist(res1)

qqnorm(res1); qqline(res1)

t.test(res1)

shapiro.test(res1)

par(mfrow=c(1,1))

acf(res1,lag = 36)

Box.test(res1, lag = 8,type ="Ljung")

#----------------------------------------------------------------------

##預測

rain11=rain[,1]

testrain =ts(rain11[(9*12+1):(10*12)],start=c(2016,1),frequency=12)

testrain1=(testrain)^(1/4)

par(mfrow=c(1,1))

rain11=rain[,1]

xreg<-as.matrix(model.matrix(~month.-1))[,-1]

model.1<-arima(raints1,order=c(0,0,0),xreg=xreg)

win.graph(width=6.5,height=3,pointsize=8)

newmonth.=season(ts(rep(1,12),start=c(2016,1),freq=12))

newtrend=(1:12)*deltat(raints1)

xreg.1 <-as.matrix(model.matrix(~newmonth.))[,-1]

plot(model.1,n.ahead=12,n1=c(2007,1),xlab="Year",pch=19,ylab='(rain)^1/4',

main="Long Term Forecasts forrain",newxreg=xreg.1,type='o')

points(testrain1,col = "blue",pch=3)

#----------------------------------------------------------------------

##季節趨勢項模型&SARIMA(0,1,1)MSE比較

a<-as.vector(plot(model.1,n.ahead=12,n1=c(2007,1),xlab="Year",pch=19,ylab='(rain)^1/4',main="LongTerm Forecasts for rain",newxreg=xreg.1,type='o')$pred)

win.graph(width=6.5,height=3,pointsize=8)

plot(m1.rain,n1=c(2007,1),n.ahead=12,col='red',xlab='Year',type='o',

ylab='(rain)^1/4',

main="Long Term Forecasts forrain")

points(testrain1,col = "blue",pch=3)

b<-as.vector(plot(m1.rain,n1=c(2007,1),n.ahead=12,col='red',xlab='Year',type='o',

ylab='(rain)^1/4',main="Long TermForecasts for rain")$pred)

y=testrain1

y

sum((a-y)^2)

sum((b-y)^2)

Correlation

#-------------------------------------------------------------------

##Pm2.5和降雨量殘差分析

respm2.5=residuals(m2.pm2.5)

resrain=residuals(model)

pm2.5rain=ts.intersect(respm2.5,resrain)

plot(pm2.5rain,main="Pm2.5 and RainResidual Plot ")

#-------------------------------------------------------------------

##CCF分析

pm2.5=pm2.5rain[,1]

rain=pm2.5rain[,2]

ccf(as.vector(pm2.5),as.vector(rain),ylab='ccf',main='x=pm2.5,y=rain')

#-------------------------------------------------------------------

##pm2.5往後移動7步後兩模型殘差時間序列圖

resrainnew=ts(rain[1:101],start=c(2007,1),end=c(2015,5),frequency=12)

respm2.5new=ts(pm2.5[8:108],start=c(2007,8),end=c(2015,12),frequency=12)

pm2.5rainnew=ts.intersect(respm2.5new,resrainnew)

plot(pm2.5rainnew,main="經過d=7移動，pm2.5月均濃度和月均雨量的殘差時間序列圖")

#-------------------------------------------------------------------

##配適回歸模型

mode1final=lm(respm2.5new~ resrainnew)

summary(mode1final)

mode1final=lm(respm2.5new~ resrainnew)

res<-summary(mode1final)$residuals

#-------------------------------------------------------------------

##殘差模型配適

plot(res, type= "o ")

acf(res)

pacf(res)

eacf(res)

#-------------------------------------------------------------------

##殘差診斷

par(mfrow=c(3,1))

plot(res,type= "o ")

hist(res)

qqnorm(res); qqline(res)

t.test(res)

shapiro.test(res)

par(mfrow=c(1,1))

acf(res,lag = 36)

Box.test(res, lag = 7,type ="Ljung")

#-------------------------------------------------------------------

##相關性圖

par(mfrow=c(1,3))

plot(as.numeric(resrain),as.numeric(respm2.5))

cor(respm2.5,resrain)

abline(lm(as.numeric(resrain)~as.numeric(respm2.5)))

plot(as.numeric(resrainnew),as.numeric(respm2.5new))

cor(respm2.5new,resrainnew)

abline(lm(as.numeric(respm2.5new)~as.numeric(resrainnew)))

plot(as.numeric(raints),as.numeric(pm2.5ts))

abline(lm(as.numeric(pm2.5ts)~as.numeric(raints)))

cor(as.numeric(pm2.5ts),as.numeric(raints))

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏3 回帖

谢谢分享

a8b3c5a 发表于 2017-3-19 17:04
pm2.5 時間序列分析 R code 參考 (成果在附件中)

Pm2.5模型建立#------------------------------------- ...

谢谢你