0
私は毎年および毎週の季節性を持つ観光地で毎日の排水の流れをモデリングしています。私は降水量、排水の接続数、回帰項などのフーリエ項を含む2つの季節パターンをモデル化し、xregを使ってauto.arimaにこれらの情報をすべて入力しました。フィット感はとても良いです。適合値と観測値との間の相関は0.986である。私が抱えている問題は、予測値が歴史的価値と一致していないということです。パターンは正しいが、下にシフトされる。どんな助けでも大歓迎です。グラフ、コード、およびデータサンプルが含まれています。前もって感謝します！私が間違ってスタックをしている場合は申し訳ありません！auto.arimaからの予測は過去のデータと一致しません
MainData <- read.delim("MainData.txt")
PrecipR <- read.delim("PrecipR.txt")
library("dplyr")
library("forecast")
library("tidyr")
library('ggplot2')
#turn rain data into a daily rate (sum of rain divided by number of hours with rain)
duplicats<-duplicated(PrecipR[,1:4])
index<-which(duplicats==FALSE)
PrecipR<-PrecipR[index,]
dsum<-aggregate(Precip~Year + Month+ Day, sum, data=PrecipR)
dcount<-aggregate(Precip~Year + Month+ Day, FUN=function(Precip){sum(Precip>0)}, data=PrecipR)
prate<-dsum$Precip/dcount$Precip
precip<-cbind(dsum, dcount[,4])
colnames(precip)[5]<-"hcount"
#combine data
dat<-MainData[order(MainData$Year, MainData$Month, MainData$Day),]
dat<-merge(dat, precip, all.x=TRUE)
dat<-dat[order(dat$Year, dat$Month, dat$Day),]
#move to beginning of cycle and attach data
dat2006<-dat[which(dat$Year==2006),]
start<-which(dat$Influent.Flow==min(dat2006$Influent.Flow))
dat<-dat[-seq(1:start),]
attach(dat)
#add holiday dummy
holiday<-rep(0, nrow(dat))
holiday[which(dat$Month == 7 & dat$Day==4)]<-1
holiday[which(dat$Month == 8 & dat$Day==18 & dat$Year==2016)]<-1
holiday[which(dat$Month == 8 & dat$Day==20 & dat$Year==2015)]<-1
holiday[which(dat$Month == 8 & dat$Day==21 & dat$Year==2014)]<-1
holiday[which(dat$Month == 8 & dat$Day==22 & dat$Year==2013)]<-1
holiday[which(dat$Month == 8 & dat$Day==16 & dat$Year==2012)]<-1
holiday[which(dat$Month == 8 & dat$Day==18 & dat$Year==2011)]<-1
holiday[which(dat$Month == 8 & dat$Day==19 & dat$Year==2010)]<-1
holiday[which(dat$Month == 8 & dat$Day==20 & dat$Year==2009)]<-1
holiday[which(dat$Month == 8 & dat$Day==21 & dat$Year==2008)]<-1
holiday[which(dat$Month == 8 & dat$Day==16 & dat$Year==2007)]<-1
holiday[which(dat$Month == 8 & dat$Day==17 & dat$Year==2006)]<-1
index<-which(holiday==1)
holiday[index+1]<-1
holiday[index-1]<-1
#model auto.arima with weekly and yearly seasons using msts and Fourier terms
flow.msts<-msts(Influent.Flow, seasonal.periods=c(7, 365.25))
y <- msts(Influent.Flow, seasonal.periods=c(7,365.25))
z <- fourier(y, K=c(2,2))
covariates.msts = cbind(Precip,
Connections,
holiday)
#main model
fit <- auto.arima(y, xreg=cbind(z,covariates.msts), seasonal=FALSE)
cor(Influent.Flow, as.numeric(fit$fitted), use = "pairwise.complete.obs")
x<-seq(1:length(Influent.Flow))
plot(x, Influent.Flow, ylim=c(0,800000), main="Daily Waste Water Flow, Observed and Modeled Values", ylab="Flow (MGD)", xlab="Daily Values, 2006 to February 2017")
lines(x, fit$fitted, col="red")
#forcast between 2/17/17 and 5/23/17
p<-precip$Precip[which(precip$Year==2017,precip$Month>=2)]
covariates.forecast = cbind(Precip<-p[seq(1:95)],
Connections<-rep(max(Connections),95),
holiday<-rep(0,95))
colnames(covariates.forecast)<-colnames(covariates.msts)
zf <- fourier(y, K=c(2,2), h=95)
fc <- forecast(fit, xreg=cbind(zf,covariates.forecast), h=95)
> head(MainData)
Date Month Day Year Influent.Flow Connections
1 1/1/06 1 1 2006 141166 484672
2 1/2/06 1 2 2006 99883 484672
3 1/3/06 1 3 2006 108132 484672
4 1/4/06 1 4 2006 88356 484672
5 1/5/06 1 5 2006 89273 484672
6 1/6/06 1 6 2006 71614 484672
> head(PrecipR)
Year Month Day Hour Precip
1 2006 1 1 0 0.03
3 2006 1 1 1 0.02
4 2006 1 1 2 0.01
5 2006 1 1 3 0.00
6 2006 1 1 4 0.00
8 2006 1 1 5 0.02
