“R Notebook3”

library(readxl)
library(xts)

Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(zoo)
library(tseries)

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

    ‘tseries’ version: 0.10-55

    ‘tseries’ is a package for time series analysis and computational finance.

    See ‘library(help="tseries")’ for details.

library(forecast)
library(astsa)

Attaching package: ‘astsa’

The following object is masked from ‘package:forecast’:

    gas

library(lmtest)
library(quantmod)

Loading required package: TTR

library(rcompanion)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

Attaching package: ‘rcompanion’

The following object is masked from ‘package:forecast’:

    accuracy

#igae <- read_excel("igae.xlsx")
# I convert the igae data frame into a time series R dataset.
# This is convenient when working with ARIMA-SARIMA models:
igae<-ts(coredata(igae$igae),start=c(1993,1),frequency=12)
# The ts function transform the dataset from a data frame to a ts object
# I can see the content of IGAE in a table format:
igae

           Jan       Feb       Mar       Apr       May       Jun       Jul       Aug       Sep       Oct
1993  55.43474  56.45697  58.90055  57.13584  57.89185  57.47547  57.90238  57.12393  58.48205  57.57943
1994  57.73260  57.41307  59.74193  60.61363  60.68776  60.46936  58.96492  60.07126  61.56862  60.75661
1995  60.28145  55.69095  57.81615  53.94180  55.49763  55.57742  54.65890  55.82789  56.93347  55.48935
1996  59.12146  58.00265  59.42397  58.28434  60.55323  59.59499  60.15472  59.57412  60.48985  61.87606
1997  61.98950  60.62307  61.06119  64.52385  64.66886  64.25285  64.76396  63.98879  66.17105  67.09512
1998  66.70703  66.23632  69.44681  67.71840  68.82410  68.48282  69.00469  67.75045  69.80193  68.82344
1999  68.30702  67.57427  71.23873  69.06278  70.35437  70.28152  70.43913  69.79525  72.35957  70.44012
2000  71.98514  72.25337  74.65102  71.37476  75.42078  74.94976  74.03275  74.84266  75.71660  74.41419
2001  73.60377  71.05126  74.40388  71.92299  74.85536  73.78798  72.88881  74.03873  73.81657  73.41285
2002  71.11273  68.96655  70.53595  74.85804  74.71278  72.86338  73.61961  73.89173  73.95660  74.95119
2003  72.57080  70.94125  73.52230  73.55614  74.75987  74.30872  74.20733  72.61894  74.30248  74.77580
2004  73.61018  72.55727  77.76094  76.01498  76.99510  77.96902  75.76327  75.86297  76.74987  76.54240
2005  74.59920  73.87880  76.06342  79.21566  79.08821  78.08797  76.10833  78.49594  79.09756  78.96851
2006  79.26468  76.80669  82.03555  79.66690  84.41093  83.27966  80.64096  82.41044  82.16254  83.52620
2007  80.89832  78.34934  83.83029  81.96804  85.63411  84.98797  83.28157  84.27601  83.26693  86.59398
2008  82.47762  81.21682  80.73307  86.80412  85.47530  85.47721  85.63097  83.08253  84.12038  86.63165
2009  76.01928  73.65056  78.41848  76.25685  76.80358  79.08988  80.18104  78.07159  80.06389  82.34449
2010  77.81779  77.04208  83.84759  82.21790  82.73190  84.09684  83.78024  82.85990  83.53129  84.51650
2011  81.05954  80.11561  87.17367  82.95717  86.03863  86.87516  86.13094  87.28997  86.72954  87.54727
2012  85.35761  85.29676  90.12011  86.67087  90.12598  89.78837  89.91496  89.44841  87.31406  91.36422
2013  87.66204  85.13913  87.34883  90.45804  91.66892  89.26224  91.10692  90.19141  88.21855  92.82476
2014  88.05532  86.92439  91.51308  90.60682  93.77173  92.24084  94.15109  91.26925  91.30502  95.61883
2015  90.61837  89.26774  94.33154  92.92938  94.01845  95.59873  96.24419  93.90688  95.48785  97.29180
2016  91.73521  92.72532  94.12730  95.74951  96.00010  97.38334  95.39440  96.88204  95.73717  97.66266
2017  95.14731  93.32814  99.88218  94.25092  99.16053  99.89463  96.51398  99.22838  95.75754  99.45651
2018  97.40601  95.31908  98.98407  98.91516 102.01301 101.24582 100.14597 101.60593  97.99756 102.34203
2019  98.52326  95.93447  99.92811  97.67383 101.33425  99.86476 100.50840 100.13507  97.49343 100.54339
2020  98.31222  95.26592  96.93029  76.33528  76.90536  86.48487  90.45817  91.53160  92.96407  96.34888
2021  93.10421  90.85215  99.78204  96.17016  98.19363  98.41539  97.30023  96.06323  93.73433  96.28978
2022  94.83383  93.68875 101.66213  98.59094 101.68172 100.26133  99.73609 101.65689  98.76533 100.62790
2023  99.27366  96.75914 104.09257 100.43890 106.18635 104.44241 102.84111 105.38761 102.31953 104.86606
2024 101.17846 101.08249 102.71041                                                                      
           Nov       Dec
1993  57.75727  60.32139
1994  61.71356  62.33790
1995  57.58134  60.10253
1996  62.31630  63.34562
1997  66.53394  68.13269
1998  68.30501  70.04768
1999  71.26505  72.41569
2000  74.31652  73.07913
2001  73.04805  72.28831
2002  73.35960  73.87354
2003  72.96502  76.15752
2004  77.59718  78.61375
2005  80.25075  81.09041
2006  82.72740  82.52557
2007  84.78503  83.55080
2008  82.90397  83.15577
2009  82.19643  83.29176
2010  86.26824  86.36191
2011  90.78646  89.20682
2012  93.41055  90.41816
2013  93.04824  92.09495
2014  94.90172  95.49843
2015  96.71978  97.50800
2016 100.84951  99.78559
2017 102.79904 100.82449
2018 104.00535 100.02001
2019 102.89807  99.93652
2020  99.29231  97.96995
2021 101.12694  99.70351
2022 105.62897 103.72796
2023 107.81851 104.75121
2024                    

plot(igae)

logigae = log(igae)
plot(logigae)

seasonaldiflog = diff(logigae,lag=12)

plot(seasonaldiflog)

adf.test(seasonaldiflog,k=1)

Warning in adf.test(seasonaldiflog, k = 1) :
  p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  seasonaldiflog
Dickey-Fuller = -5.9337, Lag order = 1, p-value = 0.01
alternative hypothesis: stationary

library(astsa)
acf2(seasonaldiflog,max.lag = 24)

     [,1] [,2] [,3]  [,4] [,5]  [,6]  [,7] [,8]  [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
ACF  0.77 0.63 0.52  0.40 0.34  0.28  0.19 0.14  0.07 -0.04 -0.12 -0.21 -0.16 -0.09 -0.06 -0.02  0.01 -0.03
PACF 0.77 0.08 0.04 -0.09 0.07 -0.01 -0.10 0.02 -0.07 -0.16 -0.08 -0.12  0.28  0.10  0.02  0.00  0.04 -0.13
     [,19] [,20] [,21] [,22] [,23] [,24]
ACF  -0.05 -0.02 -0.05 -0.06 -0.04  -0.1
PACF -0.09  0.09 -0.10 -0.14  0.04  -0.2

library(forecast)
model1 <- auto.arima(igae,
              lambda = 0, # this means that the model will transform the variable with its natural logarithm  
              d=0,D=1, # we are modeling the annual % growth month by month (not the monthly %growth)
              max.p=2, max.q=1, # setting the maximum values for p and q
              max.P = 1,max.Q = 1 # setting the max. values for P and Q
              )
coeftest(model1)

z test of coefficients:

         Estimate  Std. Error  z value  Pr(>|z|)    
ar1    1.32915524  0.20762481   6.4017 1.536e-10 ***
ar2   -0.36023358  0.18844718  -1.9116  0.055929 .  
ma1   -0.62697511  0.18350579  -3.4167  0.000634 ***
sma1  -0.82712226  0.03368356 -24.5557 < 2.2e-16 ***
drift  0.00157625  0.00025156   6.2658 3.709e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

forecast_igae <- forecast(model1, h = 21)
autoplot(forecast_igae)

#CHALLENGE 1

#Serie de ruido blanco
set.seed(123)
ruido_blanco1 <- rnorm(model1$residuals, mean = 0, sd = 1)

# Graficar la serie
plot.ts(ruido_blanco, main = "Serie de Ruido Blanco")

# Análisis de la serie
# Media y varianza
media <- mean(ruido_blanco1)
varianza <- var(ruido_blanco1)

# Autocorrelación
acf(ruido_blanco1, main = "Función de Autocorrelación del Ruido Blanco")

# Mostrar resultados
cat("Media:", media, "\n")

Media: 0.03437179 

cat("Varianza:", varianza, "\n")

Varianza: 0.9424598 

#En la gráfica del serie de tiempo se puede ver como se presenta un patrón facilmente de ver, el cual los puntos que estan distribuidos están alrededor de 0 y se puede ver como la serie si se comporta como ruido blanco. En la función de autocorrelación, muestra valores de autocorrelación de 0, significando que no hya una correlación entre los valores en diferentes tiempos.La media es de 0.03437179 siendo cercana a 0, siendo correcto con las características del ruido blanco, y la varianza da 0.94, y con una desviación estándar de 0.970824 siendo razonablemente cercana a 1.

plotNormalHistogram(model1$residuals, # series
                    prob = TRUE, # change between density and frequency
                    main = "Histograma Ruido Blanco serie de 100 obs." )

#Estos resultados muestran como la serie se comporta de manera consistente con las características del ruido blanco.

#NUEVO SIN RNOM

#Serie de ruido blanco
ruido_blanco2 <- model1$residuals

# Graficar la serie
plot.ts(ruido_blanco2, main = "Serie de Ruido Blanco")

# Análisis de la serie
# Media y varianza
media <- mean(ruido_blanco2)
varianza <- var(ruido_blanco2)
desv_st <- sqrt(varianza)

# Autocorrelación
acf(ruido_blanco2, main = "Función de Autocorrelación del Ruido Blanco")

pacf(ruido_blanco2, main = "Función de Autocorrelación del Ruido Blanco")

# Mostrar resultados
cat("Media:", media, "\n")

Media: 0.0003704806 

cat("Varianza:", varianza, "\n")

Varianza: 0.0005588996 

cat("Desv Estand:",desv_st, "\n")

Desv Estand: 0.02364106 

plotNormalHistogram(model1$residuals, # series
                    prob = TRUE, # change between density and frequency
                    main = "Histograma Ruido Blanco serie de 100 obs." )

#ALTERNATIVA 2 (UNA PRUEBA ESTADISTICA // Ljung-Box)

residuos <- residuals(model1)
plot(residuos, main = "Residuos del modelo ARIMA", ylab= "Residuos")

#Realizar la prueba de Ljung-Box con el número de lags calculado
#Ho: es que los residuos no tienen autocorrelación
ljung_box_test <- Box.test(residuos, lag=2, type= "Ljung-Box")
#Imrpimir el resultado
print(ljung_box_test)

    Box-Ljung test

data:  residuos
X-squared = 0.19247, df = 2, p-value = 0.9083

ljung_box_test <- Box.test(residuos, lag=12, type= "Ljung-Box")
print(ljung_box_test)

    Box-Ljung test

data:  residuos
X-squared = 29.576, df = 12, p-value = 0.003233

#Ruido blanco que no este correlacionado por eso se hacen las pruebas para verificar esto

#La prueba de Ljung-box verifica a independiencia de los residuos en una seride de tiempo, esta prueba evalúa si hay autocorrelación en los reisduos múltiples retornos (lags). Si los residos son ruido blanco, no deberían mostrar autocorrrelación significativa. Para la hípostesis nula los residuos no están correlacionados (con ruido blanco) y el hipotesis alternativa los residuos estan correlacionados, es decir, no son ruido blanco. En conclusión, no teneos suficiente evidencia para rechazar la hipotesis nula, los reisudos no muestran una autocorrelación significativa, lo cual es consistente con el comportamiento de ruido blanco.

#El modelo ARIM no ha capturado adecuadamente que se presento una disminución, pero que no se ha ido.

#En conclusión, no se ha superado el covid.

#CHALLENGE 2 #Interpret the calibrated ARIMA-SARIMA model. Be sure to interpret the sign and magnitude of the coefficients #Los coeficientes AR (AR1=1.32915524, AR2=-0.36023358), el AR1 tiene un coeficiente psitivo y significativo que indica una fuerte correlación positiva con el valor del período anterior, y en el AR2 tiene un coeficiente negativo y de igual manera significativo, tiene una correlación negatuva con el valor de dos períodos anteriores. En el coeficiente MA1 (-0.62697511), sugiere un coeficiente negativo, aqui incia un error positivo en el período que tiende a disminuir el valor actual. En el coeficiente MA Estacional, tiene un valor negativo (-0.82712226) y es altamente significativo, y con una deriva de (0.00157625) siendo positivo y sugiere una tendencia ascendente ligera en la serie temporal. En el gráfico se puede observar como en la línea negra representa los valofres históricos del índice de actividad económica (IGAE), en la línea azul es el pronóstico central basado en el modelo ARIMA-SARIMA, y lo sombreado representa los intervalos de predicción, mostrando una incertidumbre alrededor del pronóstico central. En este modelo se muestra como el IGAE continuará con su tendencia general ascendente, aunque contendr algunas fluctuaciones, sin embargo, muestra como la incertidumbre aumenta con el tiempo.