Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- f1 = rep(NA, 56)
- for(k in 11:56) try({
- tt=ts(t[1:k-1],start=2000,freq=4)
- sttt=decompose(ts(dat[1:k-1],start=2000,freq=4))$seasonal
- datt=dat[1:k-1]-sttt
- model=lm(datt~1+I(tt)+I(tt^2))
- trend=as.numeric(model$coefficients[1]+model$coefficients[2]*t[k]
- +model$coefficients[3]*t[k]^2)
- #sez=as.numeric(forecast(sttt,h=1)$mean)
- resid=model$res
- ar1=arima(resid,order=c(1,0,0))
- pak=as.numeric(forecast(ar1,1)$mean)
- f=trend+pak
- f1[k]=f
- }, silent = T)
- f1=ts(f1,start=2000,freq=4)
- mae=c()
- mape=c()
- for(i in 11:56){
- mae[i-10]=mean(abs(dat[11:i]-f1[11:i]),na.rm=T)
- mape[i-10]=mean(abs((dat[11:i]-f1[11:i])/dat[11:i]*100),na.rm=T)
- }
- mae=ts(na.omit(mae),start=2002.5,freq=4)
- mape=ts(na.omit(mape),start=2002.5,freq=4)
- opar=par(mfrow=c(1,2))
- plot(mae,main="MAE")
- plot(mape, main="MAPE")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement