Introducción

Este análisis se enfoca en estudiar la evolución del precio de la acción de Cemargos, una empresa destacada en la producción y comercialización de cemento en Colombia y otros países de América. Se utilizarán herramientas de series de tiempo financieras para transformar los datos históricos en un formato adecuado para el análisis. El objetivo es identificar patrones temporales, evaluar la volatilidad de la acción y realizar pronósticos utilizando modelos ARIMA y GARCH, los cuales permiten modelar tanto la tendencia como la variabilidad en los retornos.

library(dynlm)
library(ggplot2)
library(rugarch)
library(readxl)
library(xts)
library(quantmod)
library(forecast)
library(tseries)
library(FinTS)

cermargos <- read.csv("D:/Users/Daniel/Desktop/Universidad/6to semestre/series de tiempo/Cemargos Stock Price History.csv")

cermargos$Date <- as.Date(cermargos$Date, format = "%m/%d/%Y")
cermargos$Price <- as.numeric(gsub(",", "", cermargos$Price))

1. Rendimiento Diario de la Acción

Este código convierte la serie de precios de CEMARGOS en un objeto de clase xts (eXtensible Time Series), que es una estructura especializada para manejar series temporales en R. La función xts() toma los precios como valores y las fechas como índice temporal, lo que permite realizar análisis y visualizaciones propios de datos cronológicos. La función class() se usa para confirmar que la transformación se realizó correctamente.

Serie temporal indexada

Este bloque de código tiene como objetivo transformar la serie de precios históricos del CEMARGOS en una estructura adecuada para el análisis de series temporales, calcular el rendimiento diario y visualizar ambos elementos gráficamente

cermargos_xts = xts(cermargos$Price, order.by = cermargos$Date)

plot(cermargos_xts, main = "Precio Historico")

Rendimiento diario de la acción

Se calcula el rendimiento diario de la acción, es decir, la variación porcentual diaria del precio.

rend_daily = dailyReturn(cermargos_xts)

plot(rend_daily, main= "Comportamiento del Precio")

2. Investigar efectos ARCH

Modelado ARIMA del Rendimiento Diario de la Acción de CEMARGOS

Se ajustó un modelo ARIMA(0,0,2) a los rendimientos diarios, lo que indica que el comportamiento actual depende de los errores de los dos días anteriores (modelo MA(2)).

Los coeficientes de medias móviles son pequeños (ma1 = 0.0363, ma2 = 0.0537), lo que sugiere que el impacto de los días anteriores es moderado.

La varianza del error es baja (σ² = 0.0004639), lo que indica que el modelo explica bien los cambios en los rendimientos. Además, los valores negativos de AIC y BIC muestran un buen ajuste estadístico.

Por último, los residuos no presentan autocorrelación (ACF1 ≈ 0), lo que significa que el modelo captura adecuadamente la estructura de la serie sin dejar patrones sin explicar.

modelo_auto_arima = auto.arima(rend_daily)

summary(modelo_auto_arima)
## Series: rend_daily 
## ARIMA(0,0,2) with zero mean 
## 
## Coefficients:
##          ma1     ma2
##       0.0363  0.0537
## s.e.  0.0193  0.0189
## 
## sigma^2 = 0.0004639:  log likelihood = 6471.76
## AIC=-12937.53   AICc=-12937.52   BIC=-12919.85
## 
## Training set error measures:
##                        ME       RMSE       MAE MPE MAPE      MASE          ACF1
## Training set 0.0002195302 0.02153037 0.0135604 NaN  Inf 0.6913934 -0.0003252309

El modelo ARIMA(0,0,2) para rend_diario sugiere que la serie puede predecirse eficazmente considerando los errores de previsión de los dos periodos anteriores.

  • La ausencia de términos AR y de diferenciación implica que la serie es estacionaria y que su valor actual no depende directamente de sus propios valores observados en el pasado, sino más bien de perturbaciones pasadas imprevisibles.

  • Los coeficientes ma1 y ma2 significativos indican la importancia de estos errores pasados en la modelización de la serie.

H0 = No es ESTACIONARIA si el p-value es mayor que 0.05

H1 = Si es ESTACIONARIA si el p-value es menor que 0.05

  • Estacionaria p-value = 0.01
ndiffs(rend_daily)
## [1] 0
adf.test(rend_daily)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  rend_daily
## Dickey-Fuller = -13.96, Lag order = 13, p-value = 0.01
## alternative hypothesis: stationary
autoplot(modelo_auto_arima)

Supuestos sobre residuos

Hipótesis nula: No hay autocorrelación en los residuos (es decir, son ruido blanco).

p-valor = 0.03588:

El modelo ARIMA(0,0,2) ajusta bastante bien los datos.

  • Los residuos muestran leve autocorrelación

  • Como es menor a 0.05, rechazamos la hipótesis nula, lo que indica que aún hay algo de autocorrelación en los residuos, aunque leve.

checkresiduals(modelo_auto_arima) 

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,2) with zero mean
## Q* = 16.49, df = 8, p-value = 0.03588
## 
## Model df: 2.   Total lags used: 10

3. Determinar si hay Heterocedasticidad

Residuales al cuadrado

  • Dado que el p-value es muchísimo menor que 0.05, rechazamos H₀. Esto significa que hay evidencia estadísticamente significativa de autocorrelación en los errores al cuadrado.

  • Los residuos presentan heterocedasticidad condicional, lo cual indica que un modelo como ARCH o GARCH sería apropiado para modelar la volatilidad de la serie

errs_cuadra = resid(modelo_auto_arima)^2
plot(errs_cuadra, main = "Errores_cuadrado")

Box.test(errs_cuadra, lag = 5, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  errs_cuadra
## X-squared = 322.22, df = 5, p-value < 2.2e-16

Regresión con los residuales al cuadrado rezagados

  • El modelo ajusta los errores al cuadrado en función de su valor rezagado

  • El coeficiente de L(errs_cuadra, 1) es positivo y altamente significativo (1.79e-01, p < 2e-16), lo que indica que hay persistencia en la varianza: si en el periodo anterior hubo un error grande, es probable que el siguiente también lo sea.

  • Hay evidencia estadística clara de heterocedasticidad en los residuos. Esto sugiere que deberías considerar un modelo como ARCH o GARCH para capturar correctamente la dinámica de la volatilidad.

Re_residuos = dynlm(errs_cuadra ~ L(errs_cuadra,1))

summary(Re_residuos)
## 
## Time series regression with "ts" data:
## Start = 2, End = 2675
## 
## Call:
## dynlm(formula = errs_cuadra ~ L(errs_cuadra, 1))
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.015450 -0.000381 -0.000331 -0.000121  0.092606 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       3.808e-04  5.079e-05   7.496  8.9e-14 ***
## L(errs_cuadra, 1) 1.790e-01  1.903e-02   9.404  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002587 on 2672 degrees of freedom
## Multiple R-squared:  0.03204,    Adjusted R-squared:  0.03168 
## F-statistic: 88.44 on 1 and 2672 DF,  p-value: < 2.2e-16

Autocorrelación y autocorrelación parcial

  • Se observan valores significativamente altos en los primeros rezagos (especialmente en los primeros 20), ya que superan las bandas azules (límites de confianza).

  • Esto sugiere que hay dependencia en la varianza: los errores grandes tienden a seguir siendo grandes, y los pequeños, pequeños.

  • En otras palabras, hay evidencia de heterocedasticidad, lo que indica que un modelo como GARCH podría ser más adecuado para capturar esa estructura de volatilidad.

autoplot(acf(errs_cuadra, lag.max = 2434, ylim = c(-0.5,1))) +
  labs(title = "Autocorrelación parcial de los errores al cuadrado") +
  xlab("Rezagos") +
  ylab("Autocorrelación parcial")

Test de Comportamiento ARCH

Dado que el p-value es muchísimo menor a 0.05, rechazamos con fuerza la hipótesis nula.

  • Hay heterocedasticidad en la serie

  • La varianza de los errores no es constante en el tiempo

  • Hay evidencia clara de efectos ARCH -> es decir, la varianza de los errores en un momento depende de errores anteriores

arch_model = ArchTest(rend_daily, lags = 1, demean = T) # lags son los rezagos
arch_model
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  rend_daily
## Chi-squared = 87.369, df = 1, p-value < 2.2e-16

Test de Comportamiento ARCH

Dado que el p-value = 0.0001113 < 0.05, rechazamos la hipótesis nula.

Hay evidencia estadísticamente significativa de heterocedasticidad en los residuos. Específicamente, los residuos presentan efectos ARCH, lo cual significa que:

  • La varianza de los errores no es constante a lo largo del tiempo.

  • Hay periodos de alta y baja volatilidad agrupados,

arch_model = ArchTest(rend_daily, lags = 1, demean = T) # lags son los rezagos
arch_model
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  rend_daily
## Chi-squared = 87.369, df = 1, p-value < 2.2e-16

4. Modelo GARCH

ugarch1 <- ugarchspec(mean.model = list(armaOrder = c(0,0)))

ugarch1
## 
## *---------------------------------*
## *       GARCH Model Spec          *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## ------------------------------------
## GARCH Model      : sGARCH(1,1)
## Variance Targeting   : FALSE 
## 
## Conditional Mean Dynamics
## ------------------------------------
## Mean Model       : ARFIMA(0,0,0)
## Include Mean     : TRUE 
## GARCH-in-Mean        : FALSE 
## 
## Conditional Distribution
## ------------------------------------
## Distribution :  norm 
## Includes Skew    :  FALSE 
## Includes Shape   :  FALSE 
## Includes Lambda  :  FALSE
ugfit1 <- ugarchfit(spec = ugarch1, data = rend_daily)
ugfit1
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(0,0,0)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000201    0.000240  0.83735  0.40240
## omega   0.000008    0.000015  0.50688  0.61224
## alpha1  0.138252    0.011443 12.08222  0.00000
## beta1   0.858339    0.055724 15.40336  0.00000
## 
## Robust Standard Errors:
##         Estimate  Std. Error  t value Pr(>|t|)
## mu      0.000201    0.003237 0.062206  0.95040
## omega   0.000008    0.000310 0.025284  0.97983
## alpha1  0.138252    0.084833 1.629687  0.10317
## beta1   0.858339    1.143606 0.750555  0.45292
## 
## LogLikelihood : 6980.797 
## 
## Information Criteria
## ------------------------------------
##                     
## Akaike       -5.2163
## Bayes        -5.2075
## Shibata      -5.2163
## Hannan-Quinn -5.2131
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                     0.2207  0.6385
## Lag[2*(p+q)+(p+q)-1][2]    0.6042  0.6466
## Lag[4*(p+q)+(p+q)-1][5]    1.0515  0.8479
## d.o.f=0
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                  0.0001227  0.9912
## Lag[2*(p+q)+(p+q)-1][5] 2.0086158  0.6164
## Lag[4*(p+q)+(p+q)-1][9] 3.3999567  0.6915
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]    0.8061 0.500 2.000  0.3693
## ARCH Lag[5]    3.2164 1.440 1.667  0.2599
## ARCH Lag[7]    3.6443 2.315 1.543  0.4008
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.3832
## Individual Statistics:             
## mu     0.2576
## omega  0.3064
## alpha1 0.3086
## beta1  0.3677
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.07 1.24 1.6
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           0.1262 0.8996    
## Negative Sign Bias  0.6606 0.5089    
## Positive Sign Bias  0.4870 0.6263    
## Joint Effect        0.9083 0.8234    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     369.9    7.734e-67
## 2    30     473.4    8.214e-82
## 3    40     597.0   5.502e-101
## 4    50     777.1   1.103e-131
## 
## 
## Elapsed time : 0.513133
# ver coeficientes
ugfit1@fit$coef
##           mu        omega       alpha1        beta1 
## 2.013768e-04 7.838201e-06 1.382517e-01 8.583390e-01
#imprimir varianza
ug_var1 = ugfit1@fit$var
autoplot(ts(ug_var1))

# residuales
ug_resid1 = (ugfit1@fit$residuals)^2

autoplot(ts(ug_resid1))

El modelo GARCH(1,1) ajustado a los rendimientos diarios captura de manera efectiva la volatilidad condicional de la serie. Los parámetros estimados son estadísticamente significativos (excepto la media), y la suma de los coeficientes α₁ + β₁ ≈ 0.94 indica una alta persistencia en la volatilidad, lo cual es característico de series financieras.

Las pruebas diagnósticas muestran que:

  • No hay autocorrelación significativa en los residuos,

  • No hay efectos ARCH residuales, y

  • Los parámetros son estables en el tiempo.

Sin embargo, el modelo presenta deficiencias en el ajuste a los valores extremos (colas de la distribución), lo cual se evidencia en el test de bondad de ajuste de Pearson, que rechaza la hipótesis de normalidad en los residuos.

5. Predicción

# pronostico 
ug_forecast1 = ugarchforecast(ugfit1, n.ahead = 30)
ug_forecast1
## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 30
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=2024-12-30]:
##         Series   Sigma
## T+1  0.0002014 0.02155
## T+2  0.0002014 0.02169
## T+3  0.0002014 0.02184
## T+4  0.0002014 0.02198
## T+5  0.0002014 0.02212
## T+6  0.0002014 0.02226
## T+7  0.0002014 0.02239
## T+8  0.0002014 0.02253
## T+9  0.0002014 0.02267
## T+10 0.0002014 0.02280
## T+11 0.0002014 0.02293
## T+12 0.0002014 0.02306
## T+13 0.0002014 0.02319
## T+14 0.0002014 0.02332
## T+15 0.0002014 0.02345
## T+16 0.0002014 0.02358
## T+17 0.0002014 0.02370
## T+18 0.0002014 0.02383
## T+19 0.0002014 0.02395
## T+20 0.0002014 0.02407
## T+21 0.0002014 0.02420
## T+22 0.0002014 0.02432
## T+23 0.0002014 0.02444
## T+24 0.0002014 0.02455
## T+25 0.0002014 0.02467
## T+26 0.0002014 0.02479
## T+27 0.0002014 0.02490
## T+28 0.0002014 0.02502
## T+29 0.0002014 0.02513
## T+30 0.0002014 0.02524
pronostico_arima <- forecast(modelo_auto_arima, h = 30)
plot(pronostico_arima)

Modelo GARCH: Espera que la acción de CEMARGOS tenga, en promedio, una pequeña pérdida diaria constante de aproximadamente -0.012%.

Modelo ARIMA (el gráfico): Espera que, después de unos pocos días, el rendimiento diario sea prácticamente cero (ni ganancia ni pérdida promedio).

Ambos modelos sugieren que no habrá grandes cambios explosivos en el rendimiento diario; se espera que sea muy cercano a cero. El GARCH es un poco más pesimista.

Modelo GARCH: Predice que la volatilidad de CEMARGOS comenzará en un nivel (1.7%) y luego irá aumentando gradualmente durante los próximos 30 días (hasta un 2.03%). Esto significa que se espera que la acción se vuelva un poco más “nerviosa” o impredecible.

Modelo ARIMA: No nos dice cómo cambiará la volatilidad. Simplemente asume que la volatilidad se mantendrá en su nivel promedio histórico.

Ninguno de los dos modelos predice grandes saltos o caídas dramáticas en el precio para los próximos 30 días. El GARCH sugiere una leve tendencia a la baja, mientras que el ARIMA sugiere que se quedará más o menos donde está.