knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
library(ggfortify)
## Registered S3 methods overwritten by 'ggfortify':
##   method                 from    
##   autoplot.Arima         forecast
##   autoplot.acf           forecast
##   autoplot.ar            forecast
##   autoplot.bats          forecast
##   autoplot.decomposed.ts forecast
##   autoplot.ets           forecast
##   autoplot.forecast      forecast
##   autoplot.stl           forecast
##   autoplot.ts            forecast
##   fitted.ar              forecast
##   fortify.ts             forecast
##   residuals.ar           forecast
library(data.table)
## 
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
## 
##     between, first, last
library(TSA)
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(tseries)
library(knitr)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(urca)
library(seasonal)
library(dse)
## Loading required package: tfplot
## Loading required package: tframe
## 
## Attaching package: 'dse'
## The following object is masked from 'package:TSA':
## 
##     acf
## The following objects are masked from 'package:forecast':
## 
##     forecast, is.forecast
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following objects are masked from 'package:stats':
## 
##     acf, simulate
library(qqplotr)
## 
## Attaching package: 'qqplotr'
## The following objects are masked from 'package:ggplot2':
## 
##     stat_qq_line, StatQqLine
library(normtest)

arroots <- function(object)
{
  if(!("Arima" %in% class(object)) &
     !("ar" %in% class(object)))
    stop("object must be of class Arima or ar")
  if("Arima" %in% class(object))
    parvec <- object$model$phi
  else
    parvec <- object$ar
  if(length(parvec) > 0)
  {
    last.nonzero <- max(which(abs(parvec) > 1e-08))
    if (last.nonzero > 0)
      return(structure(list(
          roots=polyroot(c(1,-parvec[1:last.nonzero])),
          type="AR"),
        class='armaroots'))
  }
  return(structure(list(roots=numeric(0), type="AR"),
    class='armaroots'))
}

# Compute MA roots
maroots <- function(object)
{
  if(!("Arima" %in% class(object)))
    stop("object must be of class Arima")
  parvec <- object$model$theta
  if(length(parvec) > 0)
  {
    last.nonzero <- max(which(abs(parvec) > 1e-08))
    if (last.nonzero > 0)
      return(structure(list(
          roots=polyroot(c(1,parvec[1:last.nonzero])),
          type="MA"),
        class='armaroots'))
  }
  return(structure(list(roots=numeric(0), type="MA"),
    class='armaroots'))
}

plot.armaroots <- function(x, xlab="Real", ylab="Imaginary",
    main=paste("Inverse roots of", x$type,
          "characteristic polynomial"),
    ...)
{
  oldpar <- par(pty='s')
  on.exit(par(oldpar))
  plot(c(-1,1), c(-1,1), xlab=xlab, ylab=ylab,
       type="n", bty="n", xaxt="n", yaxt="n", main=main, ...)
  axis(1, at=c(-1,0,1), line=0.5, tck=-0.025)
  axis(2, at=c(-1,0,1), label=c("-i","0","i"),
    line=0.5, tck=-0.025)
  circx <- seq(-1,1,l=501)
  circy <- sqrt(1-circx^2)
  lines(c(circx,circx), c(circy,-circy), col='gray')
  lines(c(-2,2), c(0,0), col='gray')
  lines(c(0,0), c(-2,2), col='gray')
  if(length(x$roots) > 0)
  {
    inside <- abs(x$roots) > 1
    points(1/x$roots[inside], pch=19, col='black')
    if(sum(!inside) > 0)
      points(1/x$roots[!inside], pch=19, col='red')
  }
}
##Tema
tema_personalizado <- function(base_size = 10
)
{
  color.background = "#FFFFFF" # Chart Background
  color.grid.major = "#D9D9D9" # Chart Gridlines
  color.axis.text = "#666666" # 
  color.axis.title = "#666666" # 
  color.title = "#666666"
  color.subtitle = "#666666"
  strip.background.color = '#9999CC'
  
  ret <-
    theme_bw(base_size=base_size) +
    
    # Set the entire chart region to a light gray color
    theme(panel.background=element_rect(fill=color.background, color=color.background)) +
    theme(plot.background=element_rect(fill=color.background, color=color.background)) +
    theme(panel.border=element_rect(color=color.background)) +
    
    # Format the grid
    theme(panel.grid.major=element_line(color=color.grid.major,size=.55, linetype="dotted")) +
    theme(panel.grid.minor=element_line(color=color.grid.major,size=.55, linetype="dotted")) +
    theme(axis.ticks=element_blank()) +
    
    # Format the legend, but hide by default
    theme(legend.position="none") +
    theme(legend.background = element_rect(fill=color.background)) +
    theme(legend.text = element_text(size=base_size-3,color=color.axis.title)) +
    
    theme(strip.text.x = element_text(size=base_size,color=color.background)) +
    theme(strip.text.y = element_text(size=base_size,color=color.background)) +
    #theme(strip.background = element_rect(fill=strip.background.color, linetype="blank")) +
    theme(strip.background = element_rect(fill = "grey70", colour = NA)) +
    # theme(panel.border= element_rect(fill = NA, colour = "grey70", size = rel(1)))+
    # Set title and axis labels, and format these and tick marks
    theme(plot.title=element_text(color=color.title, 
                                  size=20, 
                                  vjust=1.25, 
                                  hjust = 0.5
    )) +
    
    theme(plot.subtitle=element_text(color=color.subtitle, size=base_size+2,  hjust = 0.5))  +
    
    theme(axis.text.x=element_text(size=base_size,color=color.axis.text)) +
    theme(axis.text.y=element_text(size=base_size,color=color.axis.text)) +
    theme(text=element_text(size=base_size, color=color.axis.text)) +
    
    theme(axis.title.x=element_text(size=base_size+2,color=color.axis.title, vjust=0)) +
    theme(axis.title.y=element_text(size=base_size+2,color=color.axis.title, vjust=1.25)) +
    theme(plot.caption=element_text(size=base_size-2,color=color.axis.title, vjust=1.25)) +
    
    # Legend  
    theme(legend.text=element_text(size=base_size,color=color.axis.text)) +
    theme(legend.title=element_text(size=base_size,color=color.axis.text)) +
    theme(legend.key=element_rect(colour = color.background, fill = color.background)) +
    theme(legend.position="bottom", 
          legend.box = "horizontal", 
          legend.title = element_blank(),
          legend.key.width = unit(.75, "cm"),
          legend.key.height = unit(.75, "cm"),
          legend.spacing.x = unit(.25, 'cm'),
          legend.spacing.y = unit(.25, 'cm'),
          legend.margin = margin(t=0, r=0, b=0, l=0, unit="cm")) +
    
    # Plot margins
    theme(plot.margin = unit(c(.5, .5, .5, .5), "cm"))
  
  ret
}

ipc <- read_excel("ipc.xlsx")
ipc<-filter(ipc, Año>2009)

Metodología Box - Jenkins

La metodología Box Jenkins consta de cuatro pasos: identificación, estimación, verificación del diagnóstico y la predicción. La serie de tiempo escogida para realizar el pronóstico es el Índice de Precios al Consumidor de Colombia desde el año 2010 hasta el año 2019. El IPC mide la evolución del precio promedio de un conjunto de bienes representativa para el consumidor, la variación porcentual entre un periodo de tiempo del índice es la inflación. Como se observa en la gráfico de la serie temporal la base del IPC es el año 2018.

ipc_ts<-ts(ipc$IPC, start = c(2010,1), frequency = 12)
autoplot(ipc_ts, ts.colour = "red", ts.linetype = 1, size = 0.1 )+
  tema_personalizado() + 
  scale_x_date(NULL,breaks = scales::breaks_width("1 years"), labels = scales::label_date("%Y"))+
    theme(axis.text.x = element_text(angle =45, vjust=0.5, size = 12), panel.grid.minor = element_blank())+
    labs(title="Evolución del IPC", 
       subtitle="2010-2020 ", 
       caption="Fuente: Departamento Nacional de Estadística DANE", 
       y="IPC", 
       x= "Año",
       color=NULL)

1. Identifiación

En el gráfico de la evolución del IPC en Colombia, se observa que la serie de tiempo tiene una tendencia positiva. Los valores que toma el IPC a través del tiempo se encuentran alrededor de la tendencia, es decir, corresponde a un proceso estocástico en tendencia. De acuerdo con Gujarati (2010), la mayor parte de las series de tiempo de datos económicos tienen orden de integración uno, es decir, necesitan diferenciarse una sola vez para convertirse en estacionarias. A simple vista es fácil concluir que la media no es estable a lo largo del tiempo, de inmediato, se puede afirmar que la media no es estacionaria. Dos de los supuestos para el pronostico en series temporales es la estacionariedad en media y la estacionariedad en varianza, en otras palabras, la media y la varianza deben de ser constantes en el tiempo.

El gráfico de la función de autocorrelación (FAC) muestra que las covarianzas de los k rezagos en el tiempo tienen un decaimiento lento hacia cero y la función de autocorrelación parcial (FACP) muestra que solo el primer rezago es significativo, indicando que la serie no es estacionaria.

ggAcf(ipc_ts, lag=20)+
  tema_personalizado() + 
  scale_x_continuous(breaks = c(1:20))+
    labs(title="Función de Autocorrelación", 
       subtitle="serie IPC", 
       y="Autocorrelaciones",
       x= "Rezago",
       color=NULL)

ipc_acf<-data.frame(("AC"),transpose(data.frame(Acf(ipc_ts, lag=10, plot="false")$acf[-1])))
names(ipc_acf)<-c("Rezago",1:10)
knitr::kable(ipc_acf,caption = "Diez primeros rezagos de la FAC",digits = 2)
Diez primeros rezagos de la FAC
Rezago 1 2 3 4 5 6 7 8 9 10
AC 0.98 0.96 0.94 0.92 0.9 0.88 0.85 0.83 0.81 0.79 
ggPacf(ipc_ts, lag=20)+
  tema_personalizado() + 
  scale_x_continuous(breaks = c(1:20))+
    labs(title="Función de Autocorrelación Parcial", 
       subtitle="Serie IPC", 
       y="Autocorrelaciones parciales", 
       x= "Rezago",
       color=NULL)

ipc_pacf<-data.frame("ACP",transpose(data.frame(Acf(ipc_ts, lag=10, type="partial", plot="false")$acf[])))
names(ipc_pacf)<-c("Rezago", 1:10)
knitr::kable(ipc_pacf,caption = "Diez primeros rezagos de la FACP",digits = 2)
Diez primeros rezagos de la FACP
Rezago 1 2 3 4 5 6 7 8 9 10
ACP 0.98 -0.02 -0.03 -0.03 -0.03 -0.02 0 0 -0.01 -0.02

Para corregir la no estacionariedad en media de la serie de tiempo se aplica diferenciación al proceso hasta que se logra la estacionariedad. Aunque, el análisis gráfico ha permitido identificar la inestabilidad de la media y la carianza a lo largo de la serie, no se tiene suficiente información para aseverar que después de aplicar diferenciación al proceso la varianza no sigue dependiendo del tiempo. Según Gujarati (2010), la aplicación de una transformación logarítmica puede solucionar los problemas de estacionariedad en varianza y se aplica antes de realizar la diferenciación de la serie. A continuación, observe el gráfico de la serie después de aplicar la transformación logarítmica y la diferenciación:

ipc_trd=diff(log(ipc_ts))
autoplot(ipc_trd, ts.colour = "skyblue", ts.linetype = 1 )+
  tema_personalizado() + 
  scale_x_date(NULL,breaks = scales::breaks_width("1 years"), labels = scales::label_date("%Y"))+
    theme(axis.text.x = element_text(angle =45, vjust=0.5, size = 12), panel.grid.minor = element_blank())+
    labs(title="Serie IPC transformada y diferenciada", 
       subtitle="", 
       y="diff(log(IPC))", 
       x= "Año",
       color=NULL)

Como se observa en la gráfica de la serie transformada y diferenciada, sus valores oscilan alrededor de una línea paralela al eje x, lo que indica que la tendencia fue eliminada.

ggAcf(ipc_trd, lag=30) +
  tema_personalizado() + 
  scale_x_continuous(breaks = seq(5,30,5))+
    labs(title="Función de Autocorrelación", 
       subtitle="Serie IPC en primeras diferencias transformadas", 
       y="Autocorrelaciones",  
       x= "Rezago",
       color=NULL)

ipc_trd_acf<-data.frame(("FAC"),transpose(data.frame(Acf(ipc_trd, lag=10, plot="false")$acf[-1])))
names(ipc_trd_acf)<-c("Rezago",1:10)
knitr::kable(ipc_trd_acf,caption = "Diez primeros rezagos de la FAC",digits = 2)
Diez primeros rezagos de la FAC
Rezago 1 2 3 4 5 6 7 8 9 10
FAC 0.65 0.28 0.02 -0.08 -0.21 -0.26 -0.19 -0.07 0.05 0.19 
ggPacf(ipc_trd, lag=30, ci=0.95)+
  tema_personalizado() + 
  scale_x_continuous(breaks = seq(5,30,5))+
    labs(title="Función de Autocorrelación Parcial", 
       subtitle="Serie IPC en primeras diferencias transformadas", 
       y="Autocorrelaciones parciales", 
       x= "Rezago",
       color=NULL)

ipc_trd_pacf<-data.frame(("FACP"),transpose(data.frame(Acf(ipc_trd, lag=10, plot="false", type="partial")$acf[])))
names(ipc_trd_pacf)<-c("Rezagos",1:10)
knitr::kable(ipc_trd_pacf,caption = "Diez primeros rezagos de la FACP",digits = 2)
Diez primeros rezagos de la FACP
Rezagos 1 2 3 4 5 6 7 8 9 10
FACP 0.65 -0.23 -0.1 0.02 -0.24 -0.02 0.07 0.01 0.05 0.19

Como se puede observar en el gráfico de la serie transformada y la función de autocorrelación, al parecer existe un componente estacional que aún no ha sido identificado. A continuación, se muestra el gráfico de la subserie estacional, en el que cada mes es separado y muestra la evolución anual de la serie permitiendo identificar si la serie de tiene algún componente estacional.

ggmonthplot(ipc_trd, year.labels = TRUE, year.labels.left = TRUE) +
  tema_personalizado() + 
    labs(title="Subserie estacional", 
       subtitle="Serie IPC en primeras diferencias transformadas", 
       y="Diff(log(IPC))",  
       x= "Meses",
       color=NULL)

Las líneas azules representan la media del conjunto de años por mes. los dos primeros meses presentan una media más alta con respecto a los demás meses, lo que indica que a la serie se le debe aplicar una diferenciación estacional.

ipc_est<-diff(ipc_trd, lag=12)
autoplot(ipc_est, ts.colour = "skyblue", ts.linetype = 1 )+
  tema_personalizado() + 
  scale_x_date(NULL,breaks = scales::breaks_width("1 years"), labels = scales::label_date("%Y"))+
    theme(axis.text.x = element_text(angle =45, vjust=0.5, size = 12), panel.grid.minor = element_blank())+
    labs(title="Serie IPC desestacionalizada", 
       subtitle="", 
       y="", 
       x= "Año",
       color=NULL)

ggAcf(ipc_est, lag=50) +
  tema_personalizado() + 
  scale_x_continuous(breaks = seq(12,50,12))+
    labs(title="Función de Autocorrelación", 
       subtitle="Serie IPC desestacionalizada", 
       y="Autocorrelaciones",  
       x= "Rezago",
       color=NULL)

ggPacf(ipc_est, lag=50, ci=0.95)+
  tema_personalizado() + 
  scale_x_continuous(breaks = seq(12,50,12))+
    labs(title="Función de Autocorrelación Parcial", 
       subtitle="Serie IPC desestacionalizada", 
       y="Autocorrelaciones parciales", 
       x= "Rezago",
       color=NULL)

En ambas funciones de autocorrelación de la serie transformada en primeras diferencias, a medida que aumentan los k rezagos la covarianza tiende a cero con un comportamiento oscilante. De acuerdo Álvarez, Calvo, Torrado y Mondragón (2013), cuando un proceso estocástico discreto presenta una atenuación sinusoidal con p rezagos significativos en la FACP y q rezagos significativos en la FAC, las estructuras son AR(p) y MA(q), es decir, los procesos generadores son AR(1) y MA(2). En la parte estacional, el análisis se realiza cada s rezagos, como el s en este caso es 12, porque la desestacionalización fue mensual se buscan los rezagos significativos cada 12 rezagos, iniciando desde el 12. Por una parte, se observa en la FAC que el rezago 12 es significativo, pero el 24 no; por lo tanto presenta una estructura SMA(1). Por otra parte, la FACP presenta 2 rezagos significativos, uno en el rezago 12 y otro en el 24; lo que indica que sigue una estructura SAR(2). En conclusión el proceso generador de datos sería un \(\ SARIMA(1, 1, 2)(2,1,1)_{12}\ \) . Sin embargo, la metodología Box-Jenkins sugiere escoger entre varios PGD el que mejor se ajuste para estimar el modelo. En la FACP puede notarse que el cuarto rezago se encuentra un poco por encima de las bandas de confianza, por ende, se podría pensar en un PGD \(\ SARIMA(1, 1, 4)(2,1,1)_{12}\ \) .

sarima1<-Arima(log(ipc_ts), order= c(1,1,2), seasonal=list(order=c(2,1,1),period=12), method = "ML")
sarima1l<-"SARIMA(1,1,2)(2,1,1)12"
sarima2<-Arima(log(ipc_ts), order= c(1,1,4), seasonal=list(order=c(2,1,1),period=12), method = "ML")
sarima2l<-"SARIMA(1,1,4)(2,1,1)12"
sarima3<-Arima(log(ipc_ts), order= c(1,1,1), seasonal=list(order=c(2,1,1),period=12), method = "ML")
sarima3l<-"SARIMA(1,1,1)(2,1,1)12"
sarima4<-Arima(log(ipc_ts), order= c(1,1,1), seasonal=list(order=c(1,1,0),period=12), method = "ML")
sarima4l<-"SARIMA(1,1,1)(1,1,0)12"
sarima5<-Arima(log(ipc_ts), order= c(1,1,0), seasonal=list(order=c(2,1,0),period=12), method = "ML")
sarima5l<-"SARIMA(1,1,0)(2,1,0)12"
sarima6<-Arima(log(ipc_ts), order= c(0,1,2), seasonal=list(order=c(0,1,1),period=12), method = "ML")
sarima6l<-"SARIMA(0,1,2)(0,1,1)12"
bic1<-BIC(sarima1, sarima2, sarima3, sarima4, sarima5, sarima6)
aic1<-AIC(sarima1, sarima2, sarima3, sarima4, sarima5, sarima6)
comp_sarima<-data.frame(c(sarima1l, sarima2l, sarima3l, sarima4l, sarima5l, sarima6l), bic1[,c(1,2)], aic1[,2])
names(comp_sarima)<-c("Modelo","Estimadores", "BIC", "AIC")
kable(comp_sarima, row.names = FALSE, round=3)
Modelo Estimadores BIC AIC
SARIMA(1,1,2)(2,1,1)12 7 -1075.188 -1094.402
SARIMA(1,1,4)(2,1,1)12 9 -1071.303 -1096.007
SARIMA(1,1,1)(2,1,1)12 6 -1079.537 -1096.006
SARIMA(1,1,1)(1,1,0)12 4 -1060.977 -1071.957
SARIMA(1,1,0)(2,1,0)12 4 -1072.665 -1083.644
SARIMA(0,1,2)(0,1,1)12 4 -1086.533 -1097.513

Según Álvarez, Calvo, Torrado y Mondragón (2013), El mejor PGD ARIMA es la que tiene menor valor del criterio de información de Akaike AIC y el Bayesiano BIC. Bajo ambos criterios el mejor modelo es el \(\ SARIMA(0, 1, 2)(0,1,1)_{12}\ \).

2. Estimación

Realizado el análisis gráfico, se procede a aplicar las pruebas formales de raíz unitaria y estacionariedad. De acuerdo con Gujarati (2010), una caminata aleatoria, un proceso de raíz unitaria y un proceso no estacionario son sinónimos. Cabe notar que, si la serie de tiempo es un proceso de raíz unitaria, en un esquema autorregresivo la varianza y la covarianza no están definidas. Entonces, el valor absoluto del coeficiente de autocovarianza del esquema AR(p) debe ser menor que uno. \[\ Y_{t}=\rho\ Y_{t-1} + u_{t}\] \[\ Y_{t}-Y_{t-1}=\rho\ Y_{t-1}-Y_{t-1} + u_{t}\] \[\ \Delta\ Y_{t}=(\rho\ -1) Y_{t} + u_{t}\] \[\ \Delta\ Y_{t}=\delta\ Y_{t} + u_{t}\] Si \(\delta\ = 0\) entonces, se tiene raíz unitaria, es decir, el proceso es no estacionario.

Prueba de Dickey-Fuller

La prueba de Dicky-Fuller (DF) utiliza el estadístico \(\tau\ \) para hallar los valor críticos y probar si \(\delta = 0\).

\[Hipótesis= \left\{ \begin{array}{lcc} H_{0}: \delta = 0\ Existe\ una\ raíz\ unitaria\\ H_{1}: \delta < 0\ La\ serie\ es\ estacionaria\ \end{array} \right.\]

Esta prueba se aplica en tres diferentes formas: \(\ Y_{t}\) es una caminata aleatoria, \(\ Y_{t}\) es una caminata aleatoria con deriva y \(\ Y_{t}\) es una caminata aleatoria con deriva alrededor de una tendencia determinística. Después se evalúan las tres pruebas de hipótesis simultáneamente para conocer si la serie es una caminata aleatoria y su tipo. La prueba DF supone que el termino de error es ruido blanco, es decir, tiene media cero, varianza constante y no correlación, para evitar errores en la aplicación de la prueba dado que, no se conoce si los residuales se encuentran auto correlacionadas se usa la prueba Dickey-Fuller Aumentada (DFA). La DFA adiciona la suma de las diferencias de i rezagos suficientes, determinados empíricamente, para eliminar la autocorrelación de los residuales. A continuación, se toma la serie del IPC sin aplicar la primera diferencia transformada asumiendo que el modelo es una caminata aleatoria sin deriva y con un número de rezagos determinados empíricamente, en este caso, tres.

summary(ur.df(ipc_ts, lags=3, type="none"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68701 -0.12339  0.00696  0.11350  0.72747 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## z.lag.1      0.0013140  0.0003228   4.071 8.43e-05 ***
## z.diff.lag1  0.8481590  0.0911615   9.304 7.82e-16 ***
## z.diff.lag2 -0.2212842  0.1171631  -1.889   0.0613 .  
## z.diff.lag3 -0.0743904  0.0927770  -0.802   0.4242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1963 on 120 degrees of freedom
## Multiple R-squared:  0.7388, Adjusted R-squared:  0.7301 
## F-statistic: 84.86 on 4 and 120 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: 4.071 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Con un nivel de significancia del 10% el \(\tau\ \) calculado (4.071) es mayor que el valor del \(\tau\ \) crítico (-1.62). Por lo tanto, no se rechaza la hipótesis nula y se concluye que la serie IPC tiene por lo menos una raíz unitaria, descartando de inmediato que el proceso tenga orden de integración cero. A continuación, se aplica la prueba asumiendo que la serie es una caminata aleatoria con deriva.

summary(ur.df(ipc_ts, lags=3, type="drift"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68538 -0.11782  0.00391  0.10991  0.72675 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0391123  0.1481011   0.264   0.7922    
## z.lag.1      0.0008699  0.0017124   0.508   0.6124    
## z.diff.lag1  0.8479491  0.0915203   9.265 1.03e-15 ***
## z.diff.lag2 -0.2210999  0.1176220  -1.880   0.0626 .  
## z.diff.lag3 -0.0731651  0.0932542  -0.785   0.4343    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1971 on 119 degrees of freedom
## Multiple R-squared:  0.5063, Adjusted R-squared:  0.4897 
## F-statistic: 30.51 on 4 and 119 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: 0.508 8.2572 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

Con un nivel de significancia del 10% el \(\tau\ \) calculado (0.508) es mayor que el valor del \(\tau\ \) crítico (-2.57). Por lo tanto, no se rechaza la hipótesis nula y se concluye que la serie IPC tiene por lo menos una raíz unitaria. A continuación, se aplica la prueba asumiendo que la serie es una caminata aleatoria con deriva y tendencia.

summary(ur.df(ipc_ts, lags=3, type="trend"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65266 -0.11310 -0.00151  0.11648  0.68714 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.685850   0.744604   2.264   0.0254 *  
## z.lag.1     -0.023386   0.010887  -2.148   0.0338 *  
## tt           0.007222   0.003202   2.255   0.0260 *  
## z.diff.lag1  0.833094   0.090229   9.233 1.31e-15 ***
## z.diff.lag2 -0.211659   0.115729  -1.829   0.0699 .  
## z.diff.lag3 -0.053440   0.092110  -0.580   0.5629    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1938 on 118 degrees of freedom
## Multiple R-squared:  0.5267, Adjusted R-squared:  0.5067 
## F-statistic: 26.26 on 5 and 118 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -2.1481 7.3889 2.6762 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

Con un nivel de significancia del 10% el \(\tau\ \) calculado (-2.1481) es mayor que el valor del \(\tau\ \) crítico (-3.13). Por lo tanto, no se rechaza la hipótesis nula concluyendo que la serie IPC tiene por lo menos una raíz unitaria y es un proceso no estacionario.

Aplicación de la prueba a la serie transformada, diferenciada y desestacionalizada.

A esta serie se le aplica un test DFA con deriva, pues, su componente de tendencia ya fue eliminado con la primera diferenciación.

summary(ur.df(ipc_est, lags=3, type="drift"))
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.009434 -0.001522  0.000290  0.001327  0.007646 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.0000612  0.0002226  -0.275  0.78395   
## z.lag.1     -0.3609232  0.1096236  -3.292  0.00135 **
## z.diff.lag1 -0.0690380  0.1136576  -0.607  0.54487   
## z.diff.lag2 -0.0673864  0.1078648  -0.625  0.53349   
## z.diff.lag3 -0.1626972  0.0966562  -1.683  0.09527 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002344 on 106 degrees of freedom
## Multiple R-squared:  0.2378, Adjusted R-squared:  0.2091 
## F-statistic: 8.269 on 4 and 106 DF,  p-value: 7.626e-06
## 
## 
## Value of test-statistic is: -3.2924 5.4428 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

Con un nivel de significancia del 5% el \(\tau\ \) calculado (-3.2924) es menor que el valor del \(\tau\ \) crítico (-2.88) se rechaza la hipótesis nula, por lo tanto, la serie IPC transformada en primeras diferencias y desestacionalizada, es un proceso estacionario.

Estimación del modelo

coeftest(sarima6)
## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ma1   0.676919   0.094762  7.1434 9.105e-13 ***
## ma2   0.330374   0.080118  4.1236 3.730e-05 ***
## sma1 -0.999604   0.164120 -6.0907 1.124e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Los estimadores del modelo \(\ SARIMA(0, 1, 2)(0,1,1)_{12}\ \) son muy significativos.

plot(maroots(sarima6),main="Raíces inversas MA")

Como se observa en el gráfico de todas las raíces inversas se encuentran dentro del círculo unitario, lo que prueba que el modelo es invertible y estacionario.

3. Validación de los supuestos del modelo

Media cero

e<-sarima6$residuals
mean(e)
## [1] -2.982359e-05
autoplot(e, colour = "seagreen")+
  tema_personalizado() + 
    labs(title="Residuales del modelo", 
       subtitle="", 
       y="",
       x= "Año",
       color=NULL)+ geom_hline(yintercept=mean(e), colour="red")

En el gráfico se observa que el término de error oscila alrededor de cero, lo que quiere decir que el término de error no afecta intercepto. La línea roja muestra la media de los residuales.

Normalidad

qmuestral <- rstandard(sarima6)
ggplot(sarima6,aes(sample=qmuestral))+
  stat_qq_band(fill="skyblue1", distribution = "norm", alpha=0.5)+
  stat_qq(col="orangered", alpha=1) + 
  stat_qq_line(col="darkblue")+
  xlab("Cuantiles Teóricos Normales") +
  ylab("Cuantiles Muestrales") + 
  labs( title ="QQ-plot",
        subtitle = "Cuantiles muestrales vs. Cuantiles Teóricos ")+
  tema_personalizado()

jarque.bera.test(residuals(sarima6))
## 
##  Jarque Bera Test
## 
## data:  residuals(sarima6)
## X-squared = 1675.3, df = 2, p-value < 2.2e-16

En el qqplot se puede notar que los valores muestrales se ajustan a la distribución teórica. Sin embargo, en la cola izquierda hay valores que se desvían notoriamente. Se realizó el test de Jarque-Bera para determinar si el modelo sigue una distribución normal, con un nivel de significancia del 5% se rechazó la hipótesis nula, por lo tanto la distribución de los residuales no se aproxima a una distribución normal.

varianza constante

ggplot(sarima6, aes(fitted.values(sarima6),e))+
  stat_smooth(method = "loess",colour="orangered2", fill="skyblue1") +
  stat_binhex(aes(alpha=..count..),bins = 70, fill="orangered2") +
  geom_hline(yintercept = 0,col ="darkblue", linetype = "dashed")+
  tema_personalizado() +
  xlab("Valores Ajustados") +
  ylab("Residuales")+
  labs( title ="Residuos vs Valores ajustados")

Al graficar los residuos frente a los valores ajustados, se puede evidenciar que los residuales se distribuyen de forma aleatoria alrededor de cero, por lo tanto no hay heterocedasticidad.

No autocorrelación

ggAcf(sarima6$residuals, lag=20)+
  tema_personalizado() + 
  scale_x_continuous(breaks = c(1:20))+
    labs(title="Función de Autocorrelación", 
       subtitle="Residuales del modelo", 
       y="Autocorrelaciones",
       x= "Rezago",
       color=NULL)

ggPacf(sarima6$residuals, lag=20)+
  tema_personalizado() + 
  scale_x_continuous(breaks = c(1:20))+
    labs(title="Función de Autocorrelación Parcial", 
       subtitle="Residuales del modelo", 
       y="Autocorrelaciones parciales",
       x= "Rezago",
       color=NULL)

Los gráficos de FAC y FACP de los residuales del modelo estimado muestran que los residuales no están autocorrelacionados.

4. Pronósitco

pronostico<-forecast::forecast(sarima6, h =5, level=95)
pronostico$mean<-exp(pronostico$mean)
pronostico$lower<-exp(pronostico$lower)
pronostico$upper<-exp(pronostico$upper)
pronostico$x<-exp(pronostico$x)
pronostico$fitted<-exp(pronostico$fitted)
pronostico$residuals<-exp(pronostico$residuals)
autoplot(pronostico, ts.colour = "red", ts.linetype = 1, size = 0.1 )+
  tema_personalizado() + 
  scale_x_date(NULL,breaks = scales::breaks_width("1 years"), labels = scales::label_date("%Y"))+
    theme(axis.text.x = element_text(angle =45, vjust=0.5, size = 12), panel.grid.minor = element_blank())+
    labs(title="Pronóstico del IPC", 
       subtitle="", 
       y="IPC", 
       x= "Año",
       color=NULL)

autoplot(ipc_ts, ts.colour = "red", ts.linetype = 1, size = 0.1 )+
  tema_personalizado() + 
  scale_x_date(NULL,breaks = scales::breaks_width("1 years"), labels = scales::label_date("%Y"))+
    theme(axis.text.x = element_text(angle =45, vjust=0.5, size = 12), panel.grid.minor = element_blank())+
    labs(title="Evolución del IPC", 
       subtitle="2010-2020 ", 
       caption="Fuente: Departamento Nacional de Estadística DANE", 
       y="IPC", 
       x= "Año",
       color=NULL)

Se puede evidenciar que la serie pronosticada presenta una trayectoria similar a la serie original. Con un nivel de significancia del 5% se espera que el IPC para el mes de septiembre se encuentre entre 104.6193 y 105.607, con un pronóstico puntual de 105.112. Se espera que haya una inflación mensual de 0.14%, una inflación año corrido de 1.26% y una inflación de 1.79% con respecto a septiembre del año anterior.

infmes<-(105.1120-104.96)*100/104.96
infanual<-(105.1120-103.26)*100/103.26
infancor<-(105.1120-103.80)*100/103.80

kable(pronostico)
Point Forecast Lo 95 Hi 95
Sep 2020 105.1120 104.6193 105.6070
Oct 2020 105.1753 104.2149 106.1446
Nov 2020 105.3083 103.9319 106.7030
Dec 2020 105.6997 104.0025 107.4246
Jan 2021 106.4494 104.4743 108.4620

5. Bibliografía

Álvarez, R. A. R., Calvo, J. A. P., Torrado, C. A. M., y Mondragón, J. A. U. (2013). Fundamentos de econometría intermedia: Teoría y aplicaciones. Bogotá D.C., Colombia: Universidad de los Andes.

Gujarati, D. N., y Porter, D. C. (2010). Econometría (5a. ed). México, D.F: McGraw-Hill Interamericana.