Advertisement
Guest User

Untitled

a guest
Oct 31st, 2014
134
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
R 0.80 KB | None | 0 0
  1. f1 = rep(NA, 56)
  2. for(k in 11:56) try({
  3.  
  4. tt=ts(t[1:k-1],start=2000,freq=4)
  5. sttt=decompose(ts(dat[1:k-1],start=2000,freq=4))$seasonal
  6. datt=dat[1:k-1]-sttt
  7.  
  8. model=lm(datt~1+I(tt)+I(tt^2))
  9.  
  10. trend=as.numeric(model$coefficients[1]+model$coefficients[2]*t[k]
  11. +model$coefficients[3]*t[k]^2)
  12.  
  13. #sez=as.numeric(forecast(sttt,h=1)$mean)
  14.  
  15. resid=model$res
  16. ar1=arima(resid,order=c(1,0,0))
  17. pak=as.numeric(forecast(ar1,1)$mean)
  18.  
  19. f=trend+pak
  20.  
  21. f1[k]=f
  22. }, silent = T)
  23.  
  24.  
  25.  
  26. f1=ts(f1,start=2000,freq=4)
  27. mae=c()
  28. mape=c()
  29.  
  30. for(i in 11:56){
  31. mae[i-10]=mean(abs(dat[11:i]-f1[11:i]),na.rm=T)
  32. mape[i-10]=mean(abs((dat[11:i]-f1[11:i])/dat[11:i]*100),na.rm=T)
  33. }
  34.  
  35. mae=ts(na.omit(mae),start=2002.5,freq=4)
  36. mape=ts(na.omit(mape),start=2002.5,freq=4)
  37.  
  38. opar=par(mfrow=c(1,2))
  39.  
  40. plot(mae,main="MAE")
  41. plot(mape, main="MAPE")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement