library(readxl)
library(forecast)
serie.pib <-
read_excel("C:/users/liizm/documents/econometria/PIB trimestral.xlsx",
col_types = c("skip", "numeric"),
skip = 5)
serie.pib.ts <- ts(data = serie.pib,
start = c(2009, 1),
frequency = 4)
serie.pib.ts %>%
autoplot(main = "PIB, El Salvador 2009-2021",
xlab = "Años/Meses",
ylab = "Indice")
ma2_4 <- ma(serie.pib.ts, 4, centre = TRUE)
autoplot(serie.pib.ts,main = "PIB, El Salvador 2009-2021",
xlab = "Años/Meses",
ylab = "Indice")+
autolayer(ma2_4,series = "Tt")
library(magrittr)
Yt <- serie.pib.ts
Tt <- ma2_4
SI <- Yt - Tt
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE)
#Promediando los resultados de cada trimestre
St <- St - sum(St) / 4
#Generar la serie de factores para cada valor de la serie original
St <-
rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 4)
autoplot(St,
main = "Factores Estacionales",
xlab = "Años/Meses",
ylab = "Factor Estacional")
Cálculo del Componente Irregular It
It<-Yt-Tt-St
autoplot(It,
main = "Componente Irregular",
xlab = "Aos/Meses",
ylab = "It")
Descomposición Aditiva (usando la libreria stats)
descomposicion_aditiva<-decompose(serie.pib.ts,type = "additive")
autoplot(descomposicion_aditiva,main="Descomposición Aditiva",xlab="Años/Meses")
Descomposición Aditiva usando libreria feasts
library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, PIB")+xlab("Años/Meses")
Componente Tendencia Ciclo [Tt=TCt]
Tt<- ma(serie.pib.ts, 4, centre = TRUE)
autoplot(Tt,main = "Componente Tendencia [Ciclo]", xlab = "Años/Meses",ylab = "Tt")
Cálculo de Factores Estacionales [St]
SI<-Yt/Tt #Serie sin tendencia.
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE)
St <- St*4/sum(St)
#Generar la serie de factores para cada valor de la serie original
St <-
rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 4)
autoplot(St,
main = "Factores Estacionales",
xlab = "Años/Meses",
ylab = "Factor Estacional")
Cálculo del Componente Irregular [It]
It<-Yt/(Tt*St)
autoplot(It,
main = "Componente Irregular",
xlab = "Años/Meses",
ylab = "It")
Descomposición Multiplicativa (usando la libreria stats)
descomposicion_multiplicatica<-decompose(serie.pib.ts,type = "multiplicative")
autoplot(descomposicion_multiplicatica,main="Descomposición Multiplicativa",xlab="Años/Meses")
Descomposición Multiplicativa usando libreria feasts
library(tsibble)
library(feasts)
library(ggplot2)
Yt %>% as_tsibble() %>%
model(classical_decomposition(value, type = "multiplicative")) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Multiplicativa, PIB") + xlab("Años/Meses")
Descomposición usando la libreria TSstudio
library(TSstudio)
ts_decompose(Yt, type = "additive", showline = TRUE)
Estimar el modelo HW
library(forecast)
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.8742302
## beta : 0
## gamma: 0
##
## Coefficients:
## [,1]
## a 7429.8128051
## b 48.3963750
## s1 0.9751646
## s2 1.0058150
## s3 0.9783521
## s4 1.0406683
#Generar el pronóstico:
PronosticosHW<-forecast(object = ModeloHW,h = 4,level = c(0.95))
PronosticosHW
## Point Forecast Lo 95 Hi 95
## 2022 Q1 7292.485 6847.413 7737.557
## 2022 Q2 7570.373 6930.977 8209.768
## 2022 Q3 7411.019 6646.235 8175.804
## 2022 Q4 7933.429 7119.875 8746.982
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()
Usando Forecast (Aproximación por Espacios de los Estados ETS )
library(forecast)
#Generar el pronóstico:
PronosticosHW2<-hw(y = Yt,
h = 4,
level = c(0.95),
seasonal = "multiplicative",
initial = "optimal")
PronosticosHW2
## Point Forecast Lo 95 Hi 95
## 2022 Q1 7284.036 6676.405 7891.668
## 2022 Q2 7541.782 6703.701 8379.863
## 2022 Q3 7515.116 6517.016 8513.215
## 2022 Q4 8063.144 6844.733 9281.555
library(TSstudio)
library(forecast)
ts_plot(Yt,Xtitle = "Años/Meses")
Orden de Integración
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 Integración") %>% kable_material()
Ordenes de Integración
|
|
x
|
|
d
|
1
|
|
D
|
1
|
#Gráfico de la serie diferenciada
Yt %>%
diff(lag = 4,diffences=D) %>%
diff(diffences=d) %>%
ts_plot(title = "Yt estacionaria")
Verificar los valores para (p,q) & (P,Q)
Yt %>%
diff(lag = 4,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36 )
Estimación del modelo SARIMA
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)[4]
##
## Coefficients:
## sar1 sma1
## -0.3980 -0.4550
## s.e. 0.3732 0.3688
##
## sigma^2 = 77377: log likelihood = -331.67
## AIC=669.34 AICc=669.9 BIC=674.89
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 6.530071 258.7683 120.9372 -0.004468918 2.054645 0.3680668
## ACF1
## Training set -0.09806694
modelo_estimado %>% autoplot(type="both")+theme_solarized()
modelo_estimado %>% check_res(lag.max = 36)
Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)
Verificación de sobre ajuste/sub ajuste
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
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)[4]> <ARIMA(0,1,0)(0,1,1)[4]> <ARIMA(0,1,0)(1,1,0)[4]>
## # ... 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 x 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_010_011 Orders 75930. -332. 668. 669. 672. <cpl [0]> <cpl [4]>
## 2 arima_010_110 Orders 78594. -333. 669. 670. 673. <cpl [4]> <cpl [0]>
## 3 arima_original Orders 77377. -332. 669. 670. 675. <cpl [4]> <cpl [4]>
## 4 arima_automatico Orders 62738. -334. 676. 677. 683. <cpl [1]> <cpl [4]>