R Tutorial

• Load the packages

library(tsibble) # convert time-series data in a tsibble
library(fable) # main package for modeling and forecasting with time-series data
library(feasts) # time-series data visualization
library(dplyr) # data manipulation

• Then load the data

data = read.csv("portal_timeseries.csv")
data_ts <- data |>
  mutate(month = yearmonth(date)) |>
  as_tsibble(index = month)
data_ts

ARIMA models

• Our simple autoregressive models were a good start
• But both the variables we looked at showed strong seasonal components that they didn’t capture
• To add seasonal components to our model we need to learn about a more complex version of AR models
• ARIMA models
• ARIMA stands for AutoRegressive Integrated Moving Average
• We’ve already learned the AutoRegressive part, so let’s talk about the other two

Moving average models

• Note that this is not using a moving average like you did when decomposing time-series
• MA models are like AR models in that they use past observations to predict the next step
• But instead of using the values themselves they use the past errors
• So a 1st order MA model looks like this

\(y_t = c + \theta_1 e_{t-1} + \epsilon_t\)

• So if \(\theta_1\) > 0 then if the observation has greater than the mean at the previous time step it is like to be greater than the mean at the current time step
• Positive MA components say that if the previous observation was an outlier then the next observation is expected to be an outlier in the same direction
• This makes sense for NDVI
• Greenness isn’t really driven by greenness in the same way a populations abundance is
• But if it’s greener than normal this month it’s probably a good year so it’s likely to be greener than normal next month

• Both AR and MA structures can be present in time-series
• We can combine them in an ARMA model

\(y_t = c + b_1 y_{t-1} + \theta_1 \epsilon_{t-1} + \epsilon_t\)

Accounting for trends

• An assumption when fitting models of this form is that the time-series is “stationary”
• Basically - there is no general trend in the data
• This is handled using “differencing” which is similar to what we did when decomposing time-series
• To take out the trend before model fitting the data is differenced
• Typically by subtracting the previous value

\(y_t' = y_t - y_{t-1}\)

• So, if we difference the data we are now modeling the change from time-step to time-step, not the actual values

• Adding AR, MA, and differencing together gives us an ARIMA model
  • AR: Autoregressive - Uses past values to predict future values
  • I: Integrated - Uses differencing to handle trends
  • MA: Moving Average - Uses past errors to predict future errors

Seasonal models

• One of the things ARIMA models can do easily in R is fit seasonal models
• The ARIMA model structure is available in fable in the ARIMA() function

arima_model = model(data_ts, ARIMA(NDVI))

• If we don’t provide it any details on model structure it will try to find the best fitting ARIMA model

report(arima_model)

• In ARIMA notation we have 3 numbers, given in parentheses, labeled (p, d, q)
• p is the AR order, d is the degree of differencing, and q is the MA order
• So this model is a 3rd order MA model, with no differencing, and no autoregressive terms

\(y_t = c + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + \theta_3 \epsilon_{t-3} + \epsilon_t\)

• This is also shown in the Coefficients table by the ma1, ma2, and ma3 terms
• So, if NDVI was above average over the last 3 months we expect it to be above average now
• This lines up with our idea that MA might be better for NDVI since it indicates a good year

• But what is the sar1 coefficient?
• That’s our season signal
• A seasonal autoregressive term of order 1
• The seasonal component is also shown in PDQ notation in the second set of parentheses
• We model season components with lags of the length of a full seasonal cycle
• In this case a full cycle is 1 year
• So the SAR1 part of the model is

\(y_t = c + b_{12} y_{t-12}\)

• The value in the current month is related to the value observed 1 year ago
• This makes a lot of sense for NDVI
• Ecosystems are typically greener during the summer than in the winter
• So a good piece of information is how green the ecosystem was at this time last year
• Our full model would therefore be

\(y_t = c + \theta_1 e_{t-1} + \theta_2 e_{t-2} + \theta_3 e_{t-3} + b_{12} y_{t-12} + \epsilon_t\)

• What does this look like?

arima_model_aug = augment(arima_model)
autoplot(arima_model_aug, NDVI) + autolayer(arima_model_aug, .fitted, color = "orange")

• Looks good
• Check the residuals

gg_tsresiduals(arima_model)

• Seasonal autocorrelation is now gone

You do: * Fit an ARIMA model to the rain data * Plot your data with the model fit on top * Plot the residuals

• What do you think about this result?