```{r} # Chapter 11 Code

set.seed(12345)

stationary_ts <- arima.sim(model = list(order = c(1,0,0), ar = 0.5 ), n = 500)

library(TSstudio)

ts_plot(stationary_ts, title = “Stationary Time Series”, Ytitle = “Value”, Xtitle = “Index”)

set.seed(12345)

non_stationary_ts <- arima.sim(model = list(order = c(1,1,0),ar = 0.3 ),n = 500)

ts_plot(non_stationary_ts, title = “Non-Stationary Time Series”, Ytitle = “Value”, Xtitle = “Index”)

data(AirPassengers)

ts_plot(AirPassengers, title = “Monthly Airline Passenger Numbers 1949-1960”, Ytitle = “Thousands of Passengers”, Xtitle = “Year”)

ts_plot(diff(AirPassengers, lag = 1), title = “AirPassengers Series - First Differencing”, Xtitle = “Year”, Ytitle = “Differencing of Thousands of Passengers”)

ts_plot(diff(diff(AirPassengers, lag = 1), 12), title = “AirPassengers Series - First and Seasonal Differencing”, Xtitle = “Year”, Ytitle = “Differencing of Thousands of Passengers”)

ts_plot(diff(log(AirPassengers), lag = 1), title = “AirPassengers Series - First Differencing with Log Transformation”, Xtitle = “Year”, Ytitle = “Differencing/Log of Thousands of Passengers”)

set.seed(12345)

library(plotly)

p1 <- plot_ly() p2 <- plot_ly()

for(i in 1:20){ rm <- NULL rw <- arima.sim(model = list(order = c(0, 1, 0)), n = 500) p1 <- p1 %>% add_lines(x = time(rw), y = as.numeric(rw)) p2 <- p2 %>% add_lines(x = time(diff(rw)), y = as.numeric(diff(rw))) }

p1 %>% layout(title = “Simulate Random Walk”, yaxis = list(title = “Value”), xaxis = list(title = “Index”)) %>% hide_legend()

p2 %>% layout(title = “Simulate Random Walk with First-Order Differencing”, yaxis = list(title = “Value”), xaxis = list(title = “Index”)) %>% hide_legend()

set.seed(12345)

ar2 <- arima.sim(model = list(order = c(2,0,0), ar = c(0.9, -0.3)), n = 500)

ts_plot(ar2, title = “Simulate AR(2) Series”, Ytitle = “Value”, Xtitle = “Index”)

md_ar <- ar(ar2) md_ar

par(mfrow=c(1,2)) acf(ar2) pacf(ar2)

set.seed(12345)

ma2 <- arima.sim(model = list(order = c(0, 0, 2), ma = c(0.5, -0.3)), n = 500)

ts_plot(ma2, title = “Simulate MA(2) Series”, Ytitle = “Value”, Xtitle = “Index”)

md_ma <- arima(ma2, order = c(0,0,2), method = “ML”) md_ma - par(mfrow=c(1,2)) acf(ma2) pacf(ma2)

set.seed(12345)

arma <- arima.sim(model = list(order(1,0,2), ar = c(0.7), ma = c(0.5,-0.3)), n = 500)

ts_plot(arma, title = “Simulate ARMA(1,2) Series”, Ytitle = “Value”, Xtitle = “Index”)

arma_md <- arima(arma, order = c(1,0,2)) arma_md

par(mfrow=c(1,2)) acf(arma) pacf(arma)

arima(arma, order = c(1, 0, 1), include.mean = FALSE)

arima(arma, order = c(2, 0, 1), include.mean = FALSE)

arima(arma, order = c(1, 0, 2), include.mean = FALSE)

arima(arma, order = c(2, 0, 2), include.mean = FALSE)

library(forecast) best_arma <- arima(arma, order = c(1, 0, 2), include.mean = FALSE)

checkresiduals(best_arma)

ar_fc <- forecast(md_ar, h = 100)

plot_forecast(ar_fc, title = “Forecast AR(2) Model”, Ytitle = “Value”, Xtitle = “Year”)

data(Coffee_Prices)

robusta_price <- window(Coffee_Prices[,1], start = c(2000, 1))

ts_plot(robusta_price, title = “The Robusta Coffee Monthly Prices”, Ytitle = “Price in USD”, Xtitle = “Year”)

acf(robusta_price)

robusta_price_d1 <- diff(robusta_price)

par(mfrow=c(1,2)) acf(robusta_price_d1) pacf(robusta_price_d1)

robusta_md <- arima(robusta_price, order = c(1, 1, 0))

summary(robusta_md)

checkresiduals(robusta_md)

data(USgas) # ——– Code Chank 39 ——– ts_plot(USgas, title = “US Monthly Natural Gas consumption”, Ytitle = “Billion Cubic Feet”, Xtitle = “Year”)

USgas_split <- ts_split(USgas, sample.out = 12)

train <- USgas_split\(train test <- USgas_split\)test

par(mfrow=c(1,2)) acf(train, lag.max = 60) pacf(train, lag.max = 60)

USgas_d12 <- diff(train, 12)

ts_plot(USgas_d12, title = “US Monthly Natural Gas consumption - First Seasonal Difference”, Ytitle = “Billion Cubic Feet (First Difference)”, Xtitle = “Year”)

USgas_d12_1 <- diff(diff(USgas_d12, 1))

ts_plot(USgas_d12_1, title = “US Monthly Natural Gas consumption - First Seasonal and Non-Seasonal Differencing”, Ytitle = “Billion Cubic Feet (Difference)”, Xtitle = “Year”)

par(mfrow=c(1,2)) acf(USgas_d12_1, lag.max = 60) pacf(USgas_d12_1, lag.max = 60)

p <- q <- P <- Q <- 0:2

arima_grid <- expand.grid(p,q,P,Q) names(arima_grid) <- c(“p”, “q”, “P”, “Q”) arima_grid\(d <- 1 arima_grid\)D <- 1

arima_grid$k <- rowSums(arima_grid)

library(dplyr)

arima_grid <- arima_grid %>% filter(k <= 7) # ——– Code Chank 49 ——– arima_search <- lapply(1:nrow(arima_grid), function(i){ md <- arima(train, order = c(arima_grid\(p[i], 1, arima_grid\)q[i]), seasonal = list(order = c(arima_grid\(P[i], 1, arima_grid\)Q[i])), method = “CSS”) # Cambiar el método de optimización si es necesario results <- data.frame(p = arima_grid\(p[i], d = 1, q = arima_grid\)q[i], P = arima_grid\(P[i], D = 1, Q = arima_grid\)Q[i], AIC = md$aic) }) %>% bind_rows() %>% arrange(AIC)

head(arima_search)

USgas_best_md <- arima(train, order = c(1,1,1), seasonal = list(order = c(2,1,1))) USgas_best_md

USgas_test_fc <- forecast(USgas_best_md, h = 12)

accuracy(USgas_test_fc, test)

test_forecast(USgas, forecast.obj = USgas_test_fc, test = test)

final_md <- arima(USgas, order = c(1,1,1), seasonal = list(order = c(2,1,1)))

checkresiduals(final_md)

USgas_fc <- forecast(final_md, h = 12)

plot_forecast(USgas_fc, title = “US Natural Gas Consumption - Forecast”, Ytitle = “Billion Cubic Feet”, Xtitle = “Year”)

USgas_auto_md1 <- auto.arima(train) USgas_auto_md1

USgas_auto_md2 <- auto.arima(train, max.order = 5, D = 1, d = 1, ic = “aic”, stepwise = FALSE, approximation = FALSE) USgas_auto_md2

df <- ts_to_prophet(AirPassengers) %>% setNames(c(“date”, “y”))

df\(lag12 <- dplyr::lag(df\)y, n = 12)

library(lubridate)

df\(month <- factor(month(df\)date, label = TRUE), ordered = FALSE)

df$trend <- 1:nrow(df)

par <- ts_split(ts.obj = AirPassengers, sample.out = 12)

train <- par$train

test <- par$test

train_df <- df[1:(nrow(df) - 12), ] test_df <- df[(nrow(df) - 12 + 1):nrow(df), ]

md1 <- tslm(train ~ season + trend + lag12, data = train_df)

checkresiduals(md1)

md2 <- auto.arima(train, xreg = cbind(model.matrix(~ month,train_df)[,-1], train_df\(trend, train_df\)lag12), seasonal = TRUE, stepwise = FALSE, approximation = FALSE)

summary(md2)

checkresiduals(md2)

fc1 <- forecast(md1, newdata = test_df)

fc2 <- forecast(md2, xreg = cbind(model.matrix(~ month,test_df)[,-1], test_df\(trend, test_df\)lag12))

accuracy(fc1, test) accuracy(fc2, test) ```