# Estimating ARIMA Models using bayesforecast package # Here, I'm going to try to predict the internal temperature of a cow in the morning # Load packages #### library(bayesforecast) library(fma) library(ggplot2) # Read in data #### # This is the internal temperature of a cow in the morning (degrees F) ipc <- cowtemp # Plot the data and identify any unusual observations #### autoplot(object = ipc,main = "Daily Morning Temperature of a Cow",ylab="Temperature (degrees F)") # This time series isn't stationary, so we should look at the correlation plots. #### g1 = autoplot(object = diff(ipc),main = "Differenced Cow Temp",y = "Temperature (degrees F") g2 = ggacf(y = diff(ipc)) g3 = ggpacf(y = diff(ipc)) gridExtra::grid.arrange(g1,g2,g3, layout_matrix = matrix(c(1,2,1,3),nrow = 2)) # The variance doesn't look stable over time. Do we need a Box-Cox transformation? # What type of model should we use? #### # p is the order of the AR term - 2 based on PACF plot? # q is the order of the MA term - Maybe not needed or 1? # d is the number of differcing required to make the time series stationary - 1? # I used a seasonal ARIMA because I thought that is what we would need for the flu data # and it looked like there may be some periodicity in the cow data # An ARIMA model might have been better though sf1 = stan_sarima(ts = ipc,order = c(1,1,0),seasonal = c(1,0,1), prior_sar = beta(2,2),prior_sma = beta(2,2),chains = 3) # I fit this with beta priors because it was suggested by the documentation, # but I don't understand why those are the priors of choice when using seasonal # components. The regular ARIMA model uses normal priors, which makes sense. summary(sf1) # Fit using auto.sarima? #### sf2 <- auto.sarima(ts = ipc, chains = 3, iter = 5000, adapt.delta = 0.9) summary(sf2) # This suggests maybe a moving average model instead # Check the adjustment of each parameter #### mcmc_plot(object = sf1) # Look at model fit #### autoplot(sf1)+labs(title = "Posterior Predict", y="Temp") # It looks pretty good to me! # Check the model residuals #### check_residuals(sf1) # This looks OK, but there are still strange things happening with seasonal patterns # Maybe I need a better seasonal model, or no seasonal model at all. # Predict for the next 25 days #### autoplot(object = forecast(sf1,h = 25),ylab="Temp") # My model autoplot(object = forecast(sf2,h = 25),ylab="Temp") # auto.sarima