Sección 8

Capitulo 8

# Chapter 8 Code

# -------- Code Chank 1 --------
library(TSstudio)

## Warning: package 'TSstudio' was built under R version 4.2.3

data(USgas)

ts_info(USgas)

##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10

# Using the series time index to set the start and end point of each partiton
train <- window(USgas, 
                start = time(USgas)[1], 
                end = time(USgas)[length(USgas) - 12])

test <- window(USgas, 
               start = time(USgas)[length(USgas) - 12 + 1], 
               end = time(USgas)[length(USgas)])

ts_info(train)

##  The train series is a ts object with 1 variable and 226 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2018 10

ts_info(test)

##  The test series is a ts object with 1 variable and 12 observations
##  Frequency: 12 
##  Start time: 2018 11 
##  End time: 2019 10

# -------- Code Chank 2 --------
# The sample.out argument set the size of the testing partition
# (and therefore the training partition)
USgas_partitions <- ts_split(USgas, sample.out = 12)

train <- USgas_partitions$train
test <- USgas_partitions$test

ts_info(train)

##  The train series is a ts object with 1 variable and 226 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2018 10

ts_info(test)

##  The test series is a ts object with 1 variable and 12 observations
##  Frequency: 12 
##  Start time: 2018 11 
##  End time: 2019 10

# -------- Code Chank 3 --------
library(forecast)

## Warning: package 'forecast' was built under R version 4.2.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

md <- auto.arima(train)
# -------- Code Chank 4 --------
checkresiduals(md)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.95, df = 18, p-value = 0.1263
## 
## Model df: 6.   Total lags used: 24

# -------- Code Chank 5 --------
fc <- forecast(md, h = 12)
# -------- Code Chank 6 --------
accuracy(fc, test)

##                     ME      RMSE      MAE       MPE     MAPE      MASE
## Training set  5.844136  97.81626 73.42657 0.1170672 3.522348 0.6376860
## Test set     37.847885 103.22848 81.46603 1.3107987 3.261643 0.7075062
##                      ACF1 Theil's U
## Training set -0.004183172        NA
## Test set     -0.046708926 0.3404092

# -------- Code Chank 7 --------
test_forecast(actual = USgas,
              forecast.obj = fc,
              test = test) 

# -------- Code Chank 8 --------
library(forecast)

naive_model <- naive(train, h  = 12)
test_forecast(actual = USgas,
              forecast.obj = naive_model,
              test = test)

accuracy(naive_model, test)

##                      ME     RMSE      MAE        MPE     MAPE     MASE
## Training set  -1.028444 285.6607 228.5084 -0.9218463 10.97123 1.984522
## Test set     301.891667 499.6914 379.1417  9.6798015 13.28187 3.292723
##                   ACF1 Theil's U
## Training set 0.3761105        NA
## Test set     0.7002486  1.499679

# -------- Code Chank 9 --------
snaive_model <- snaive(train, h = 12)
test_forecast(actual = USgas,
              forecast.obj = snaive_model,
              test = test)

accuracy(snaive_model, test)

##                    ME     RMSE      MAE      MPE     MAPE     MASE       ACF1
## Training set 33.99953 148.7049 115.1453 1.379869 5.494048 1.000000  0.4859501
## Test set     96.45000 164.6967 135.8833 3.612060 5.220458 1.180103 -0.2120929
##              Theil's U
## Training set        NA
## Test set     0.4289964

# -------- Code Chank 10 --------
md_final <- auto.arima(USgas)

fc_final <- forecast(md_final, h = 12)
# -------- Code Chank 11 --------
plot_forecast(fc_final,
              title = "The US Natural Gas Consumption Forecast",
              Xtitle = "Year",
              Ytitle = "Billion Cubic Feet")

# -------- Code Chank 12 --------
fc_final2 <- forecast(md_final, 
                      h = 60, 
                      level = c(80, 90))

plot_forecast(fc_final2,
              title = "The US Natural Gas Consumption Forecast",
              Xtitle = "Year",
              Ytitle = "Billion Cubic Feet")

# -------- Code Chank 13 --------
fc_final3 <- forecast_sim(model = md_final,
                          h = 60, 
                          n = 500) 
# -------- Code Chank 14 --------
library(plotly)

## Warning: package 'plotly' was built under R version 4.2.3

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

fc_final3$plot %>% 
  layout(title = "US Natural Gas Consumption - Forecasting Simulation",
         yaxis = list(title = "Billion Cubic Feet"),
         xaxis = list(title = "Year"))

# -------- Code Chank 15 --------
set.seed(1234)

#USgas_forecast$summary_plot