###################################### # ARIMA Models: U.S. Cattle Markets # ###################################### remove(list=ls()) # Load packages library(readxl) library(fpp3) library(tsibble) library(tsibbledata) library(tidyverse) library(fable) library(dplyr) library(lubridate) library(ggplot2) library(cowplot) library(ggpubr) library(feasts) # Read data setwd("C:/Documents/Classes/Econ4115/Code") data <- read_excel(na.omit("cattle.xlsx"), sheet = 1, col_names = TRUE) head(data) tail(data) # Create a time series object data1 <- data%>% filter(YEAR > 1929) %>% mutate(YEAR == 1930:2022) %>% as_tsibble(index = YEAR) ############################# # Regular Time Series Plots # ############################# # Regular time series plot of cattle inventories data1 %>% autoplot(TOTAL_COWS) + ggtitle("Total U.S. Inventory of Cows (Jan 1)") + ylab("Thousands of Cows") + xlab("Years") # Regular time series plot of cattle prices data1 %>% autoplot(PRICE_ALL_CATTLE) + ggtitle("Average Price of All Cattle") + ylab("Dollars per CWT") + xlab("Years") ############################# # Non-Seasonal ARIMA Models # ############################# ################# # Cattle Prices # ################# # Select the best ARIMA model automatically fitprices <- data1 %>% model(ARIMA(PRICE_ALL_CATTLE ~ PDQ(0,0,0))) report(fitprices) # Forecast and plot fitprices %>% forecast(h = 10) %>% autoplot(slice(data1,(n()-40):n())) # ACF and PACF data1 %>% ACF(PRICE_ALL_CATTLE) %>% autoplot() data1 %>% ACF(difference(PRICE_ALL_CATTLE)) %>% autoplot() data1 %>% PACF(difference(PRICE_ALL_CATTLE)) %>% autoplot() # The automatic search indicates an ARIMA(1,1,2) w/ drift. # MANUAL SEARCH: # The ACF is exponentially decaying which indicates non-stationarity. # The ACF of the differenced data has significant spikes at lags 2 & 5. # The PACF also has (near) significant spikes at first 5 lags. # Cattle prices may follow an ARIMA(5,1,5) process. # Search the best non-seasonal ARIMA model fitprices1 <- data1 %>% model(ARIMA(PRICE_ALL_CATTLE ~ pdq(0:5, 1, 0:5) + PDQ(0, 0, 0), stepwise = FALSE, approximation = FALSE)) report(fitprices1) # A more thorough search indicates an ARIMA(4,1,1) w/ drift. ############################ # Total Cattle Inventories # ############################ # Select the best ARIMA model automatically fitcows <- data1 %>% model(ARIMA(TOTAL_COWS ~ PDQ(0,0,0))) report(fitcows) # Forecast and plot fitcows %>% forecast(h=10) %>% autoplot(slice(data1,(n()-90):n())) # ACF and PACF data1 %>% ACF(TOTAL_COWS) %>% autoplot() data1 %>% ACF(difference(TOTAL_COWS)) %>% autoplot() data1 %>% PACF(difference(TOTAL_COWS)) %>% autoplot() # The automatic search indicates an ARIMA(4,1,0). # MANUAL SEARCH: # The ACF is exponentially decaying which indicates non-stationarity. # The ACF of the differenced data has significant (cyclical) spikes at 5 of the first 6 lags. # The PACF has (near) significant spikes at first 3 lags. # Cattle inventories may follow an ARIMA(3,1,0) process. # Search the best non-seasonal ARIMA model fitcows1 <- data1 %>% model(ARIMA(TOTAL_COWS ~ pdq(0:5, 1, 0:5) + PDQ(0, 0, 0), stepwise = FALSE, approximation = FALSE)) report(fitcows1)