library(readxl)
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
serie.ivae <- read_excel("C:/Users/lupita nieto/Downloads/Practica .xlsx",col_types = c("skip", "numeric"),
skip = 0)
serie.ivae.ts<- ts(data = serie.ivae, start = c(2009,1),
frequency = 12)
Yt<-serie.ivae.ts
Grafico De Serie Estacionaria
library(TSstudio)
library(forecast)
ts_plot(Yt,Xtitle= "Años/Meses")
Identificacion
Orden De Integracion
library(kableExtra)
library(magrittr)
d<-ndiffs(Yt)
D<-nsdiffs(Yt)
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()
Ordenes de Integracion
|
|
x
|
|
d
|
1
|
|
D
|
1
|
# Gráfico de la serie diferenciada
Yt %>%
diff(lag=12,diffences=D) %>%
diff(diffences=d) %>%
ts_plot(title= "Yt estacionaria")
valores para (p,q) & (P,Q)
Yt %>%
diff(lag=12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36)
Estimacion Del Modelo
library(forecast)
library(ggthemes)
modelo_estimado<-Yt %>%
Arima(order = c(0, 1, 0),
seasonal = c(1, 1, 1))
summary(modelo_estimado)
## Series: .
## ARIMA(0,1,0)(1,1,1)[12]
##
## Coefficients:
## sar1 sma1
## -0.1022 -0.7232
## s.e. 0.1266 0.1135
##
## sigma^2 = 6.779: log likelihood = -351.22
## AIC=708.43 AICc=708.6 BIC=717.38
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05202814 2.477766 1.68543 0.02782791 1.663587 0.4366552
## ACF1
## Training set 0.06837002
modelo_estimado %>% autoplot(type="both")+theme_solarized()
modelo_estimado %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
library(forecast)
#Estimar el modelo
modeloHW<-HoltWinters(x = Yt,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
modeloHW
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.8472467
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 118.3719784
## b 0.1600947
## s1 0.9529095
## s2 1.0192349
## s3 1.0440936
## s4 0.9962326
## s5 1.0009929
## s6 0.9838373
## s7 0.9554392
## s8 1.0145286
## s9 1.0872315
## s10 0.9748891
## s11 0.9626701
## s12 1.0059813
#Generando el pronostico:
pronosticosMW<-forecast(object = modeloHW,h=12, level=c(0.95))
pronosticosMW
## Point Forecast Lo 95 Hi 95
## Apr 2022 112.9503 107.46492 118.4358
## May 2022 120.9752 113.57246 128.3779
## Jun 2022 124.0929 115.22236 132.9634
## Jul 2022 118.5640 108.86876 128.2592
## Aug 2022 119.2908 108.50116 130.0804
## Sep 2022 117.4038 105.81297 128.9947
## Oct 2022 114.1680 101.97016 126.3657
## Nov 2022 121.3911 107.66977 135.1125
## Dec 2022 130.2643 114.88360 145.6450
## Jan 2023 116.9603 102.34981 131.5708
## Feb 2023 115.6485 100.48405 130.8129
## Mar 2023 121.0126 83.19827 158.8270
#Gráfico comparativo
Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-pronosticosMW$fitted
grafico_comparativo<-cbind(Yt, Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)
Verificacion de sobre ajuste/sub ajuste
library(tsibble)
##
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library(feasts)
## Loading required package: fabletools
##
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
##
## accuracy, forecast
library(fable)
library(fabletools)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
a<-Yt %>% as_tsibble() %>%
model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)),
arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a)
## # A mable: 1 x 4
## arima_original arima_010_011 arima_010_110
## <model> <model> <model>
## 1 <ARIMA(0,1,0)(1,1,1)[12]> <ARIMA(0,1,0)(0,1,1)[12]> <ARIMA(0,1,0)(1,1,0)[12]>
## # … with 1 more variable: arima_automatico <model>
a %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla
tabla
## # A tibble: 4 × 9
## `Model name` .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
## 1 arima_automatico Orders 6.13 -347. 702. 702. 714. <cpl [1]> <cpl [12]>
## 2 arima_010_011 Orders 6.72 -352. 707. 707. 713. <cpl [0]> <cpl [12]>
## 3 arima_original Orders 6.78 -351. 708. 709. 717. <cpl [12]> <cpl [12]>
## 4 arima_010_110 Orders 7.99 -361. 726. 726. 731. <cpl [12]> <cpl [0]>
(Cross Validated)
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt<-Yt %>% as_tsibble() %>% rename(IVAE=value)
data.cross.validation<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
TSCV<-data.cross.validation %>%
model(ARIMA(IVAE ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>%
forecast(h=1) %>% accuracy(Yt)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 1 observation is missing at 2022 abr.
print(TSCV)
## # A tibble: 1 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA(IVAE ~ pdq(0, … Test 0.0506 2.94 2.09 0.00243 2.00 0.542 0.494 0.155