#Instalación de librerias
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr   1.1.4     ✔ readr   2.1.5
## ✔ forcats 1.0.0     ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0     ✔ tibble  3.2.1
## ✔ purrr   1.0.2     ✔ tidyr   1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(astsa)
library(foreign)
library(timsac)
library(vars)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Loading required package: sandwich
## 
## Attaching package: 'strucchange'
## 
## The following object is masked from 'package:stringr':
## 
##     boundary
## 
## Loading required package: urca
## Loading required package: lmtest
library(lmtest)
library(mFilter)
library(dynlm)
library(nlme)
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(lmtest)
library(broom)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(knitr)
library(MASS)
library(parallel)
library(car)
library(mlogit)
## Loading required package: dfidx
## 
## Attaching package: 'dfidx'
## 
## The following object is masked from 'package:MASS':
## 
##     select
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(dplyr)
library(tidyr)
library(forecast)
## 
## Attaching package: 'forecast'
## 
## The following object is masked from 'package:nlme':
## 
##     getResponse
## 
## The following object is masked from 'package:astsa':
## 
##     gas
library(fpp2)
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──
## ✔ fma       2.5     ✔ expsmooth 2.3
## ── Conflicts ───────────────────────────────────────────────── fpp2_conflicts ──
## ✖ forecast::getResponse() masks nlme::getResponse()
## ✖ car::some()             masks purrr::some()
## 
## Attaching package: 'fpp2'
## 
## The following object is masked from 'package:astsa':
## 
##     oil
library(stats)
library(quantmod)
## Loading required package: xts
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'xts'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Loading required package: TTR
library(openxlsx)
#Se quitaron las fechas del excel, ya que no se leian bien en R
#Para leer el excel
db = read.xlsx("metal.xlsx")
db
##     Metal
## 1   15136
## 2   16733
## 3   20016
## 4   17708
## 5   18019
## 6   19227
## 7   22893
## 8   23739
## 9   21133
## 10  22591
## 11  26786
## 12  29740
## 13  15028
## 14  17977
## 15  20008
## 16  21354
## 17  19498
## 18  22125
## 19  25817
## 20  28779
## 21  20960
## 22  22254
## 23  27392
## 24  29945
## 25  16933
## 26  17892
## 27  20533
## 28  23569
## 29  22417
## 30  22084
## 31  26580
## 32  27454
## 33  24081
## 34  23451
## 35  28991
## 36  31386
## 37  16896
## 38  20045
## 39  23471
## 40  21747
## 41  25621
## 42  23859
## 43  25500
## 44  30998
## 45  24475
## 46  23145
## 47  29701
## 48  34365
## 49  17556
## 50  22077
## 51  25702
## 52  22214
## 53  26886
## 54  23191
## 55  27831
## 56  35406
## 57  23195
## 58  25110
## 59  30009
## 60  36242
## 61  18450
## 62  21845
## 63  26488
## 64  22394
## 65  28057
## 66  25451
## 67  24872
## 68  33424
## 69  24052
## 70  28449
## 71  33533
## 72  37351
## 73  19969
## 74  21701
## 75  26249
## 76  24493
## 77  24603
## 78  26485
## 79  30723
## 80  34569
## 81  26689
## 82  26157
## 83  32064
## 84  38870
## 85  21337
## 86  19419
## 87  23166
## 88  28286
## 89  24570
## 90  24001
## 91  33151
## 92  24878
## 93  26804
## 94  28967
## 95  33311
## 96  40226
## 97  20504
## 98  23060
## 99  23562
## 100 27562
## 101 23940
## 102 24584
## 103 34303
## 104 25517
## 105 23494
## 106 29095
## 107 32903
## 108 34379
## 109 16991
## 110 21109
## 111 23740
## 112 25552
## 113 21752
## 114 20294
## 115 29009
## 116 25500
## 117 24166
## 118 26960
## 119 31222
## 120 38641
## 121 14672
## 122 17543
## 123 25453
## 124 32683
## 125 22449
## 126 22316
## 127 27595
## 128 25451
## 129 25421
## 130 25288
## 131 32568
## 132 35110
## 133 16052
## 134 22146
## 135 21198
## 136 19543
## 137 22084
## 138 23816
## 139 29961
## 140 26773
## 141 26635
## 142 26972
## 143 30207
## 144 38687
## 145 16974
## 146 21697
## 147 24179
## 148 23757
## 149 25013
## 150 24019
## 151 30345
## 152 24488
## 153 25156
## 154 25650
## 155 30923
## 156 37240
## 157 17466
## 158 19463
## 159 24352
## 160 26805
## 161 25236
## 162 24735
## 163 29356
## 164 31234
## 165 22724
## 166 28496
## 167 32857
## 168 37198
## 169 13652
## 170 22784
## 171 23565
## 172 26323
## 173 23779
## 174 27549
## 175 29660
## 176 23356
#Convertir los datos en una serie temporal mensual a partir de neero del 1/31/1980. Este es un paso crucial para asegurar que los datos esten en el formato adecuado para el análisis de series temporales.

metal_ts <- ts(db, start=c(1980,1), frequency = 12) #cuenta con un frecuencia mensual.
#Se empleo el autoarima paara ajustar automaticamente un modelo ARIMA a la serie temporal. 

autoarima <- auto.arima(metal_ts)
summary(autoarima)
## Series: metal_ts 
## ARIMA(1,1,2)(0,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2     sma1
##       0.4299  -1.4673  0.5339  -0.6600
## s.e.  0.2984   0.2658  0.2340   0.0799
## 
## sigma^2 = 5399312:  log likelihood = -1497.05
## AIC=3004.1   AICc=3004.48   BIC=3019.57
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE    MAPE      MASE       ACF1
## Training set -170.7415 2208.571 1619.979 -1.403046 6.55824 0.8234203 0.01025396
plot(metal_ts,  main="Serie de tiempo", ylab="Metal", col="pink")

#El modelo ajustado fue un ARIMA(1,1,2)(0,1,1)[12]
#Inicialmente se realizo la prueba de Dickey Fuller para verificar la estacionariedad de la serie de tiempo original.

adf.test(metal_ts)
## Warning in adf.test(metal_ts): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  metal_ts
## Dickey-Fuller = -6.9449, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
# El resultado muestra un valor de Dickey-Fuller =-6.9449 y un p-value de 0.01, lo que sugiere que la serie es estacionaria al nivel de 1%. 
#A pesar de que indique es es una serie estacionaria, se aplica una diferenciación a la serie para asegurarse de que sea estacionaria, aunque visulamente si parece serlo. Esto es un buen enfoque para mejorar el modelo.

metaldiff1 = diff(metal_ts) 
adf.test(metaldiff1) 
## Warning in adf.test(metaldiff1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  metaldiff1
## Dickey-Fuller = -8.9901, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
plot(metaldiff1) 

#Se confirma que la serie diferenciada es estacionaria.
#se utilizaron las funciones acf y pacf para determinar los valores apropiados de los parámetros "p" y "q" del modelo. 

acf(metaldiff1) #función de autocorrelación #q = 2 ó 3

#Existe una relación inversa ya que las lineas tambien se extienden en la parte posterior.

pacf(metaldiff1)  #función de autocorrelación parcial #p = 4 ó 5

#Entonces el (p,d,q) será (4,1,2) ó (5,1,3) 
#modelo = arima(Prices.ts,order=c(p,d,q)) 
# Al inicio trate de usar otra serie de numeros pero eran muy altos, asi que lo dividi por grupos.
#Se ajusta manualmente un modelo ARIMA (4,1,2).

modelo1 = arima(metal_ts, order=c(4,1,2))
summary(modelo1)
## 
## Call:
## arima(x = metal_ts, order = c(4, 1, 2))
## 
## Coefficients:
##          ar1      ar2      ar3     ar4      ma1      ma2
##       0.0138  -0.2173  -0.0325  0.1290  -0.8128  -0.1131
## s.e.  0.2438   0.0836   0.1008  0.0786   0.2353   0.2160
## 
## sigma^2 estimated as 24651531:  log likelihood = -1738.73,  aic = 3491.46
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE     MASE        ACF1
## Training set 591.259 4950.906 3878.712 -1.029695 15.61999 0.802427 -0.00951974
#El modelo muestra un MAPE (Error Absoluto Medio Porcentual) de 15.61999, lo cual indica un ajuste relativamente alto.
#Se hara un segundo modelo para ARIMA(5,1,3) ver si se puede reducir el margen de error.

modelo2 = arima(metal_ts, order=c(5,1,3))
## Warning in log(s2): NaNs produced
summary(modelo2)
## 
## Call:
## arima(x = metal_ts, order = c(5, 1, 3))
## 
## Coefficients:
##           ar1      ar2     ar3      ar4     ar5     ma1      ma2      ma3
##       -1.3983  -0.8497  -0.266  -0.1683  0.1523  0.7623  -0.6145  -0.8908
## s.e.   0.0859   0.1445   0.160   0.1418  0.0827  0.0465   0.0476   0.0418
## 
## sigma^2 estimated as 20488528:  log likelihood = -1724.06,  aic = 3466.12
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 551.6064 4513.548 3633.964 -0.6872361 14.69134 0.7517937
##                     ACF1
## Training set -0.02627683
#Este modelo muestra un MAPE (Error Absoluto Medio Porcentual) de 14.69134, lo cual indica un ajuste relativamente alto pero menor en comparación al primer modelo.
#Finalmente usando el modelo ARIMA (5,1,3) se realizó un pronóstico para los próximos 5 meses.

pronostico=forecast::forecast(modelo2,h=5) #los doble puntos son para usar una función en especifico de una libreria
plot(pronostico)

pronostico
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Sep 1994       27261.98 21461.03 33062.94 18390.19 36133.78
## Oct 1994       26335.24 20161.97 32508.51 16894.04 35776.44
## Nov 1994       25093.40 18800.06 31386.75 15468.57 34718.24
## Dec 1994       27960.67 21665.88 34255.46 18333.62 37587.72
## Jan 1995       23635.84 17231.04 30040.65 13840.54 33431.15
#Observamos que el pronostico muestra una ligera caida inicial en los valores y luego una pequeña recuperación, lo que sigue el patron estacional y la tendencia de la serie. Los intervalos de confianza son relativamente altos, ya que son del 80% y 95%, lo cual proporciona un rango estimado de posibles valores futuros.

Conclusiones

  1. El proceso de selección de modelos arima utilizando tanto el autoarima, como el análisis manual de ACF y PACF permitió identificar un modelo que proporciona uun ajuste preciso a los datos. El modelo final ARIMA (5,1,3) logró un MAPE más bajo, sugiriendo que este es más adecuado para la serie temporal dada.

  2. La confirmación de la estacionariedad es vital para el éxito del modelo de ARIMA, ya que estos modelos asumen que la serie temporal es estacionaria.

  3. El pronóstico de 5 meses permite realizar proyeccciones futuras con un nivel de confianza especifico, lo cual es útil para la toma de decisiones basadas en datos históricos del metal, pudiendo ver que las predicciones siguen la tendencia esperada de la serie temporal.

  4. Al comparar el auto.arima con el el modelo ajustado manualmente, podemos visualizar que el autoarima tiene un margen de error mas bajo (MAPE 6.55) que el modelo que ajuste manualmente (14.69). Esto puede ser debido a que autoarima selecciona el modelo basadonse en otros criterios como el AIC, mientras que yo uso criterios como el p-value y el MAPE. Aunque el modelo ARIMA (5,1,3) ofrece un ajuste adecuado, siempre hay espacio para mejoras adicionales para que no se sobreajuste. Recomendaciones: manteniendo el modelo lo más simple posible, usar la prueba de Dickey-Fuller para determinar si la serie ya es estacionaria antes de aplicar una diferenciación adicional y evaluar su rendimiento.

  5. El proceso de raíz visualizado en la imagen es un proceso estacionario con tendencia ascendente.