Parcial 3 - MMF
Isabella Hernandez O
1. Descripción de la base de datos
La siguiente base de datos se obtuvo de https://www.investing.com/commodities/us-soybeans-historical-data en ella, se tiene el historial mensual desde febrero de 1990 hasta noviembre del 2022 sobre el precio futuro de la soya. Estos datos se obtienen por medio de La Asociación Estadounidense de la Soya (ASA) la cual representa a los productores de soya de EE.UU en cuestiones de política nacional e internacional importantes para la industria de la soya. Esta base de datos llamada los futuros de soja (Soybeans Futures) están disponibles para negociar en la Bolsa de Comercio de Chicago (CBOT®)
library(readxl)
soya <- read_excel("C:/Users/Lenovo/Downloads/SEPTIMOSEMESTRE/MMF/Parcial3MMF/USSoyBeans.xlsx")
View(soya)2. Definir la serie de tiempo
Teniendo en cuenta que los datos utilizados estan distribuidos regularmente, es decir que tenemos un dato cada año desde febrero de 1990 hasta el mes de noviembre del 2022 entonces, en el siguiente fragmento de código definimos el objeto de la serie de tiempo ts, teniendo en cuenta que los datos son mensuales e inician en el año 1990 es decir que su frecuencia es de 12.
soya.ts <- ts(soya$Price , start = c(1990,02), frequency=12)
head(soya.ts) Feb Mar Apr May Jun Jul
1990 566.50 595.00 632.50 607.25 623.50 594.50frequency(soya.ts)[1] 12plot(soya.ts)
plot(decompose(soya.ts))
El análisis nos revela que tenemos que eliminar la estacionalidad de la serie.
3. Si es necesario, utilizar transformación Box-Cox / logaritmos para estabilizar la varianza.
No me permite realizar el BoxCox.
4. Modelar la serie de tiempo elegida como función de una tendencia, tendencia cuadrática y/o medias estacionales. Presente el diagnóstico del modelo elegido.
modelo1<-lm(soya.ts~time(soya.ts))
summary(modelo1)
Call:
lm(formula = soya.ts ~ time(soya.ts))
Residuals:
Min 1Q Median 3Q Max
-360.0 -202.2 -11.9 134.6 743.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -48329.421 2435.558 -19.84 <2e-16 ***
time(soya.ts) 24.516 1.214 20.20 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 228.4 on 392 degrees of freedom
Multiple R-squared: 0.5099, Adjusted R-squared: 0.5087
F-statistic: 407.9 on 1 and 392 DF, p-value: < 2.2e-16En esta primer regresión podemos ver que nuestro intercepto es de -48329.421 y nuestra variable independiente es de 24.516 lo que significa que cada año que pasa aumenta en 24.516 el precio de venta de la soya. Aunque las variables son significativas, note que el valor R-cuadrado es de 0.50 por lo tanto, la estimación no es muy eficiente se puede tener una mejor.
A continuación se hizo un modelo con un ajuste cuadrático para mejorar el ajuste.
modelo2<-lm(soya.ts~time(soya.ts)+ I(time(soya.ts)^2))
summary(modelo2)
Call:
lm(formula = soya.ts ~ time(soya.ts) + I(time(soya.ts)^2))
Residuals:
Min 1Q Median 3Q Max
-385.70 -189.11 -6.68 124.09 757.00
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.602e+05 5.749e+05 1.670 0.0957 .
time(soya.ts) -9.808e+02 5.731e+02 -1.711 0.0878 .
I(time(soya.ts)^2) 2.505e-01 1.428e-01 1.754 0.0802 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 227.8 on 391 degrees of freedom
Multiple R-squared: 0.5138, Adjusted R-squared: 0.5113
F-statistic: 206.6 on 2 and 391 DF, p-value: < 2.2e-16Con este ajuste vemos que el tiempo, y el tiempo al cuadrado no son significativos, ademas las variables no son significativas.
5. En caso de estacionalidad de la serie de tiempo, aplicar diferencias estacionales.
Para ello, verificamos cuantas diferencias estacionales se deben realizar.
ndiffs(soya.ts)[1] 1Como se puede ver, es necesario aplicar diferenciación una vez con el fin de lograr que la serie sea estacionaria.
dif.soya.ts <- diff(soya.ts)
plot(dif.soya.ts)
6. Usar prueba Dickey-Fuller para evaluar el orden de integración de la serie.
Anteriormente, realizamos una vez el proceso de diferenciación tal como se pedia por ende, verificaremos si finalmente la serie es ESTACIONARIA por medio de la pruebaDicker-Fuller ‘adf.test’
\(H_0\) = La serie es no estacionaria.
\(H_1\) = La serie es estacionaria
adf.test(dif.soya.ts, alternative = c("stationary", "explosive"))Warning in adf.test(dif.soya.ts, alternative = c("stationary", "explosive")) :
p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: dif.soya.ts
Dickey-Fuller = -8.5377, Lag order = 7, p-value = 0.01
alternative hypothesis: stationaryCon la prueba Dickey-Fuller comparamos el p-value, el cual es de 0.01 con \(\alpha = 0.05\), por lo tanto se rechaza la hipótesis nula de que no es estacionaria, por lo que ya podemos empezar a crear modelos.
7. Usar herramientas ACF, PACF y/o criterios de Akaike/ Bayes para construir propuestas de modelos.
par(mfrow=c(1,2))
acf(dif.soya.ts)
pacf(dif.soya.ts)
Note que algunos rezagos se salen de los limites que tienen las ACF y PACF.
8. Ajuste del ARIMA
A continuación, por medio de la función auto.arima() se realizara un ajuste automático de un modelo ARIMA a la serie de tiempo que estamos trabajando, en este caso se obtiene ARIMA(0,0,0) es decir, que solo depende del ruido blanco.
modelo<-auto.arima(dif.soya.ts)
modeloSeries: dif.soya.ts
ARIMA(0,0,0) with zero mean
sigma^2 = 4846: log likelihood = -2225.13
AIC=4452.26 AICc=4452.27 BIC=4456.24residuo <- residuals(modelo)
acf(residuo)
Se puede concluir que los residuos no tienen solo un comportamiento de ruido blanco, por lo que el modelo arima podria no estar siendo el adecuado para realizar predicciones.
9. Realizar el pronóstico ARIMA para dos años.
A continuación, se realizan las predicciones del precio futuro de la soya para los proximos 12 meses y se visualizan graficamente teniendo de teniendo en cuenta el modelo generado anteriormente.
pred<-forecast(modelo,h=12)
plot(forecast(modelo))
Note que la predicción graficamente se ve como un valor constante y esto puede estar sucediendo porque los valores estan siendo muy cercanos unos a los otros.
10. Utilizando metodos de pronósticos simples y/o suavizamiento exponencial generar un pronóstico para dos años. Elegir el método que presente la raíz del error medio al cuadrado más pequeño.
Para poder realizar el pronostico de nuestra serie, se realizo un recorte de los ultimos dos años para pronosticar nuestra serie mediante los metodos mean(media), naive y naive estacional.
DIVPRO<-window(dif.soya.ts,start=2020,end=c(2022,11))
autoplot(DIVPRO) 
plot(SOYA2)
plot(SOYA3)
plot(SOYA4)
Una vez comparado con los dos años siguientes de la serie, concluimos que el metodo de naive estacional presenta un comportamiento mas logico y razonable para los siguientes años.
11. Modelo GARCH
Dado que se tiene una serie de tiempo financiera considero que se podria hacer uso de otro modelo que permita tener un mejor pronóstico teniendo en cuenta que la información adicional de pasado podria estar afectando la variación del pronostico y no se estaria ajustando bien a los datos.
Por medio de las grafica de ACF y PACF, algunos rezagos se salen, esto implica que la varianza es heterocedastica sin embargo, esto lo identificaremos por medio de la prube ARCHTEST.
\(H_0\) = No hay efectos ARCH.
\(H_1\) = Hay efectos ARCH
require(rugarch)
require(FinTS)
ipcArchTest <- ArchTest(dif.soya.ts, lags=1, demean = TRUE)
ipcArchTest
ARCH LM-test; Null hypothesis: no ARCH effects
data: dif.soya.ts
Chi-squared = 24.536, df = 1, p-value = 7.294e-07Dado que el p-value es menor a 0.05 rechazamos \(H0\) y por ende, la varianza es heterocedastica. A continuación, se intenta ajustar un modelo GARCH al conjunto de datos con el ARIMA que se viene trabajando. De aquí se tiene que el modelo GARCH adecuado es GARCH(1,1)
ug.spec = ugarchspec(mean.model = list(armaOrder=c(0,0)))
ugfit = ugarchfit(spec = ug.spec, data = dif.soya.ts )Predicciones modelo GARCH
pred<-ugarchforecast(ugfit,n.ahead=10)
plot(pred, which=1)
ugfit
*---------------------------------*
* 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 1.15857 2.011382 0.57601 0.564611
omega 170.83594 62.509117 2.73298 0.006276
alpha1 0.34009 0.070140 4.84880 0.000001
beta1 0.65891 0.057927 11.37477 0.000000
Robust Standard Errors:
Estimate Std. Error t value Pr(>|t|)
mu 1.15857 2.088479 0.55474 0.579071
omega 170.83594 68.287771 2.50171 0.012360
alpha1 0.34009 0.068275 4.98126 0.000001
beta1 0.65891 0.055666 11.83679 0.000000
LogLikelihood : -2132.85
Information Criteria
------------------------------------
Akaike 10.875
Bayes 10.915
Shibata 10.874
Hannan-Quinn 10.891
Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
statistic p-value
Lag[1] 1.258 0.2620
Lag[2*(p+q)+(p+q)-1][2] 1.521 0.3561
Lag[4*(p+q)+(p+q)-1][5] 2.119 0.5906
d.o.f=0
H0 : No serial correlation
Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
statistic p-value
Lag[1] 0.0007906 0.9776
Lag[2*(p+q)+(p+q)-1][5] 0.4403817 0.9665
Lag[4*(p+q)+(p+q)-1][9] 1.7360791 0.9341
d.o.f=2
Weighted ARCH LM Tests
------------------------------------
Statistic Shape Scale P-Value
ARCH Lag[3] 0.02256 0.500 2.000 0.8806
ARCH Lag[5] 0.69225 1.440 1.667 0.8258
ARCH Lag[7] 1.67778 2.315 1.543 0.7851
Nyblom stability test
------------------------------------
Joint Statistic: 0.9519
Individual Statistics:
mu 0.03124
omega 0.71250
alpha1 0.24520
beta1 0.42925
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
------------------------------------
Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
group statistic p-value(g-1)
1 20 32.80 0.02534
2 30 43.87 0.03775
3 40 42.52 0.32191
4 50 54.71 0.26678
Elapsed time : 0.181515 Los coeficientes están siendo significativos por ende, el modelo GARCH(1,1) se puede utilizar y este tendrá en cuenta información del pasado que podría estar afectando la variación del pronóstico. En este caso del modelo GARCH, se tiene en cuenta que la varianza es heterocedastica y que depende de los residuales al cuadrado rezagados como también de la varianza rezagada.
Conclusion
Aunque la predicción del modelo ARIMA(0,0,0) gráficamente se ve como un valor constante, esto puede estar sucediendo porque los valores están siendo muy cercanos unos a los otros considero que el modelo no será muy eficiente puesto que la serie inicial no era estacionaria, se le realizo una diferenciación y además, no logre utilizar transformación Box-Cox / logaritmos para estabilizar la varianza lo cual puede influir en el resultado del pronóstico de un modelo ARIMA. Sin embargo, los residuales se están saliendo de los límites establecidos, hay rezagos. Por ende, recurrir a otro tipo de modelos que permitan capturar este tipo de comportamientos de la serie, es lo más adecuado con el fin de obtener un mejor pronostico.