library(readxl)
library(TSA)
library(ggplot2)
library(ggfortify)
require(tseries)
require(forecast)
require(timsac)
require(changepoint)

Contextualización

El azucar es materia prima de muchos productos industriales, entre ellos chocolates,bebidas gaseosas, jarabes, etc. Es por esto que una variacion en los precios puede afectar la cadena alimentaria de las personas, ya sea por que se encarecen los precios finales de los productos o por lo contrario, los hace más economicos para el alcance de la gente con menores ingresos.Los factores que más afectan Su precio es en gran medida el clima y las oscilacionesde la demanda a nivel mundial.

La base de datos para este documento fue obtenida en Macrotrends, y es sobre los precios del azúcar en Estados Unidos dados en días y se encuentran en un periodo de 19 años (2002-2021). La información está dada en dolares por libra de azúcar.

Estadística descriptiva

library(readxl)
Data <- read_excel("C:/Users/miche/Downloads/sugar-prices-historical-chart-data.xlsx")
Precios <- ts(Data$VALUE,start = c(2002,1), frequency = 260)

Medidas de tendencia central y de posición

summary (Precios)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0536  0.1124  0.1413  0.1497  0.1784  0.3531

Al observar el summary del precio del azúcar, se tiene que el minimo precio ha sido 0.0536 (Septiembre 2004) y su máximo 0.3531 (Febrero 2011).También se puede concluir que la media es mayor que la mediana, por lo tanto, los datos están sesgados a la derecha.

Análisis de patrones y observaciones atípicas

plot (Precios,xlab= "Periodo", ylab= "Precio")

Analizando la gráfica, se puede notar que entre 2009 y 2013 hay una observación atípica porque hubo un gran alza en los precios. En este periodo de tiempo, los precios subieron por el temor de una escasez inminente, pues se preveía que el mundo consumiría más de lo que producía. Adicionalmente se puede decir que la serie presenta estacionalidad.

fit<-decompose(Precios, type = "multiplicative")
autoplot(fit)

De acuerdo con la descomposición,la serie de tiempo no cuenta con una tendencia, pues los valores están relativamente estables.

BoxCox.ar(Precios)
## Warning in arima0(x, order = c(i, 0L, 0L), include.mean = demean): possible
## convergence problem: optim gave code = 1

Después de realizar la función BoxCox, se concluyó que como el lambda se acerca a 0, es necesario hacer una transformación logaritmica de los datos.

Log_precios <- log (Precios)

Modelar la serie de tiempo

modelo1<-lm(Log_precios~time(Log_precios))
summary(modelo1)
## 
## Call:
## lm(formula = Log_precios ~ time(Log_precios))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.72909 -0.26338 -0.04706  0.22256  0.96331 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -59.512904   1.878501  -31.68   <2e-16 ***
## time(Log_precios)   0.028612   0.000934   30.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3435 on 4785 degrees of freedom
## Multiple R-squared:  0.164,  Adjusted R-squared:  0.1638 
## F-statistic: 938.4 on 1 and 4785 DF,  p-value: < 2.2e-16

Realizando regresión lineal,se tiene que por cada día se aumenta en 0.028 el precio del azúcar, y que ambos valores son significativos. Además, su R cuadrado es del 16,4% lo que nos indica que el modelo no se ajusta lo suficiente a los datos.

modelo2<-lm(Log_precios~time(Log_precios)+ I(time(Log_precios)^2))
summary(modelo2)
## 
## Call:
## lm(formula = Log_precios ~ time(Log_precios) + I(time(Log_precios)^2))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55545 -0.17792 -0.05913  0.15735  0.71743 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -3.738e+04  5.836e+02  -64.06   <2e-16 ***
## time(Log_precios)       3.714e+01  5.803e-01   64.01   <2e-16 ***
## I(time(Log_precios)^2) -9.227e-03  1.443e-04  -63.96   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2522 on 4784 degrees of freedom
## Multiple R-squared:  0.5493, Adjusted R-squared:  0.5491 
## F-statistic:  2915 on 2 and 4784 DF,  p-value: < 2.2e-16

Para el segundo modelo, se tiene que su R cuadrado es del 54,9% lo que nos indica que el modelo actual ajusta mucho mejor la base de datos.

ggplot(Log_precios, aes(y=Precios, x=time(Log_precios)))+geom_point()+geom_smooth(se=FALSE)

Gráfica después del ajuste logaritmico.

autoplot(modelo2, which = 1:6, ncol = 2, label.size = 3)

En esta gráfica, podemos observar los residuales de nuestra serie, los cuales son ligeramente diferentes con respecto a los ajustados. En cuanto a la gráfica que indica si se ajusta a un comportamiento normal, se puede notar que lo hace parcialmente.

Prueba Shapiro-Wilk

Se realiza la prueba de Shapiro-Wilk para saber si el componente estocástico está normalmente distribuido.

shapiro.test(rstudent(modelo2))
## 
##  Shapiro-Wilk normality test
## 
## data:  rstudent(modelo2)
## W = 0.9527, p-value < 2.2e-16

Puesto a que su p-valor es muy pqueño, se rechaza la hipotesis nula, lo que indica que el componente estocástico no está normalmente distribuido.

Función de autocorrelación

ggAcf(rstudent(modelo2))

Después de realizar la función de correlación, se puede observar que este modelo presenta muchos errores estandar.

Garantizar estacionalidad en la serie

Se realiza una prueba de Dickey-Fuller para saber si la serie de tiempo es estacionaria:

adf.test(Log_precios, alternative = c("stationary", "explosive"))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Log_precios
## Dickey-Fuller = -2.3546, Lag order = 16, p-value = 0.4282
## alternative hypothesis: stationary

Analizando el valor de p (mayor a 0.05), se concluye que la serie no es estacionaria, por lo que es necesario diferenciar.A continuación se identificará cuantas veces se debe diferenciar:

ndiffs(Log_precios)
## [1] 1

El resultado es tan solo una diferenciación.

Log_precios_ts <- diff(Log_precios)

Teniendo ya la serie diferenciada, se vuelve a realizar la prueba de Dickey-Fuller:

adf.test(Log_precios_ts, alternative = c("stationary", "explosive"))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Log_precios_ts
## Dickey-Fuller = -15.692, Lag order = 16, p-value = 0.01
## alternative hypothesis: stationary

Puesto que el p-value es menor a 0.05, se puede concluir que la serie ya es estacionaria.

##Pruebas ACF y PACF

par(mfrow=c(1,2))
acf(Log_precios_ts)
pacf(Log_precios_ts)

Después de realizar el ACF, se puede decir que la mayoría de sus valores no son significativos, porque se encuentran dentro de las lineas y azules, pero hay valores que son la excepción como la tercera linea que muestra una correlación muy positiva. En cuanto al PACF, se puede mencionar que sus valores tampoco son significativos, y al igual que en el ACF, se encuentran algunas excepciones.

ARIMA

Precios_AutoArima<-auto.arima(Log_precios_ts)
Precios_AutoArima
## Series: Log_precios_ts 
## ARIMA(0,0,0) with zero mean 
## 
## sigma^2 = 0.0004021:  log likelihood = 11919.11
## AIC=-23836.22   AICc=-23836.22   BIC=-23829.75

La función auto.arima nos sugiere un modelo ARIMA (0, 0, 0).Esto puede ser interpretado como un ruido blanco,lo que significa que los errores no están correlacionados a lo largo del tiempo.

checkresiduals(Precios_AutoArima) 

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with zero mean
## Q* = 549.26, df = 520, p-value = 0.1811
## 
## Model df: 0.   Total lags used: 520

Presente la ecuación final y describa

El resultado de este caso, es el modelo ARIMA (0, 0, 0), y este, será utilizado para realizar el pronóstico ARIMA para dos años.

Pronostico1<-Arima(Log_precios,order = c(0,0,0))
Pronostico1
## Series: Log_precios 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##          mean
##       -1.9681
## s.e.   0.0054
## 
## sigma^2 = 0.1411:  log likelihood = -2104.39
## AIC=4212.77   AICc=4212.77   BIC=4225.72
autoplot(forecast(Pronostico1))

##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.

En primer lugar, para realizar el pronóstico, y ver si realmente se ajusta, se redujo los años de la serie, desde 2015, hasta 2017

PRONOSTICO<-window(Log_precios,start=2015,end=c(2017,260))
autoplot(PRONOSTICO)

En segundo lugar, se llevaran a cabo las predicciones por medio de los métodos media y naive.

PRONOSTICO_MEAN <- meanf(PRONOSTICO,h=520)
PRONOSTICO_NAIVE <- naive(PRONOSTICO,h=520)
plot(PRONOSTICO_MEAN)

plot(PRONOSTICO_NAIVE)

Conclusiones

Después de realizar todo este proceso (pronóstico y vista de residuos) se pudo ver gráficamente que ninguna de las dos predicciones, se ajusta realmente a los datos verdaderos dados en los tiempos 2018-2020. Por tanto, el modelo no capturó ni la dinámica, ni la tendencia, pues simplemente dio como resultado un solo valor para los años utilizados.

Los pronósticos de este caso no son nada confiables ni realistas porque los resultados fueron muy diferentes a la realidad. Esto pudo haberse dado gracias a que la función auto.arima, dio como resultado un ARIMA(0, 0, 0) el cual respresenta a un ruido blanco (una constante)

Finalmente, una alternativa puede ser realizar otro tipo de tranformaciones a la serie que permitan encontrar un ARIMA diferente a un ruido blanco, y evaluar en las diferentes pruebas y funciones para saber si el modelo realmente está prediciendo correctamente.