library(ggplot2)
library(forecast)
library(tseries)
library(readxl)
SNG_GDP <- read_excel("C:/Users/samorapa/OneDrive/Estudios/universidad/Economía Internacional/Proyecto de Aula 2/Anexos y Consolidados/Consolidados Fianles/data segundo proyecto de aula.xlsx")
names(SNG_GDP)
## [1] "year" "SGP"  "CHN"  "JPN"  "KOR"  "MYS"  "HKG"  "IDN"
summary(SNG_GDP$SGP)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -3.680   4.903   7.867   7.428  10.308  15.240
library(ggplot2)
plot_GDP <- ggplot(SNG_GDP, aes(y=SGP, x=year, colour="PIB de Singapur"))+geom_line()+ggtitle("Evolución del PIB de Singapure")+ylab("PIB")+xlab("Años")
plot_GDP

#Limpiando datos
library(forecast)
SNG_GDP$SGP_ts <- ts(SNG_GDP[, c("SGP")])
library(tseries)
library(forecast)
library(ggplot2)
SNG_GDP$ma5 <- ma(SNG_GDP$SGP_ts, order=5)
SNG_GDP$ma15 <- ma(SNG_GDP$SGP_ts, order = 15)


ggplot()+geom_line(data = SNG_GDP, aes(x=year, y=SGP, colour="PIB de Singapur"))+geom_line(data = SNG_GDP, aes(x=year, y=ma5, colour="Promedio Movil de 5 periodos"))+geom_line(data = SNG_GDP, aes(x=year, y=ma15, colour="Promedio Movil de 15 periodos"))+ggtitle("Tendencia del PIB de Singapure")+xlab("Año")+ylab("PIB de Singapur")
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 14 rows containing missing values (geom_path).

library(tseries)
library(forecast)
GDP_ma5 <- ts(na.omit(SNG_GDP$ma5), frequency=12)
decomp <- stl(GDP_ma5, s.window = "periodic")
deseasonal_GDP <- seasadj(decomp)
plot(decomp)

adf.test(GDP_ma5, alternative = "stationary")
## Warning in adf.test(GDP_ma5, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  GDP_ma5
## Dickey-Fuller = -6.4567, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary
Acf(GDP_ma5, main="")

Pacf(GDP_ma5, main="")

GDP_D1 <- diff(deseasonal_GDP, differences = 1)
plot(GDP_D1)

#Test de Dickey-Fuller
adf.test(GDP_D1, alternative = "stationary")
## 
##  Augmented Dickey-Fuller Test
## 
## data:  GDP_D1
## Dickey-Fuller = -3.2453, Lag order = 3, p-value = 0.0898
## alternative hypothesis: stationary
Acf(GDP_D1, main="ACF para Primeras Diferencias")

Pacf(GDP_D1, main="Pafc para Primeras Diferencias")

auto.arima(deseasonal_GDP, seasonal = FALSE)
## Series: deseasonal_GDP 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -0.6199  -0.9230  0.9031  0.8553
## s.e.   0.0722   0.0636  0.1093  0.0983
## 
## sigma^2 estimated as 0.7268:  log likelihood=-63.47
## AIC=136.95   AICc=138.28   BIC=146.61
fit <- auto.arima(deseasonal_GDP, seasonal = FALSE)
tsdisplay(residuals(fit), lag.max = 18, main = "(2,1,2) Residuos del Modelo")

fcast <- forecast(fit, h=5)
plot(fcast)

#Probando el modelo sin estacionalidad
hold <- window(ts(deseasonal_GDP), start=50)
fit_no_hold <- arima(ts(deseasonal_GDP[-c(50:55)]), order=c(2,1,2))
fcast_no_hold <- forecast(fit_no_hold, h=5)
plot(fcast_no_hold)
lines(ts(deseasonal_GDP))

fit_w_seasonality <- auto.arima(deseasonal_GDP, seasonal = TRUE)
fit_w_seasonality
## Series: deseasonal_GDP 
## ARIMA(0,1,2)(0,0,1)[12] 
## 
## Coefficients:
##          ma1      ma2     sma1
##       0.2183  -0.3669  -0.8234
## s.e.  0.1395   0.1627   0.5500
## 
## sigma^2 estimated as 0.6806:  log likelihood=-67.18
## AIC=142.36   AICc=143.23   BIC=150.09
seas_fcast <- forecast(fit_w_seasonality, h=5)
plot(seas_fcast)

tsdisplay(residuals(fit_w_seasonality), lag.max = 18, main = "(0,1,2)(0,0,1)[12] Residuos del Modelo Estacional")

#Probando el modelo estacional
hold <- window(ts(deseasonal_GDP), start=50)
fit2_no_hold <- arima(ts(deseasonal_GDP[-c(50:55)]), order=c(0,1,2),seasonal = list(order=c(0,0,1), period=12))
fcast2_no_hold <- forecast(fit2_no_hold, h=5)
plot(fcast2_no_hold)
lines(ts(deseasonal_GDP))

fit_w_seasonality2 <- arima(deseasonal_GDP,order = c(0,1,5), seasonal = list(order=c(0,0,1), period=12))
fit_w_seasonality2
## 
## Call:
## arima(x = deseasonal_GDP, order = c(0, 1, 5), seasonal = list(order = c(0, 0, 
##     1), period = 12))
## 
## Coefficients:
##          ma1     ma2     ma3      ma4      ma5     sma1
##       0.1085  0.0967  0.3114  -0.2953  -0.5991  -0.2681
## s.e.  0.1088  0.1471  0.1648   0.1550   0.1343   0.2191
## 
## sigma^2 estimated as 0.6525:  log likelihood = -64.45,  aic = 142.9
seas_fcast2 <- forecast(fit_w_seasonality2, h=5)
plot(seas_fcast2)

tsdisplay(residuals(fit_w_seasonality2), lag.max = 18, main = "(0,1,5)(0,0,1)[12] Residuos del Modelo Estacional Ajustando Promedio Móvil")

fit_w_seasonality3 <- arima(deseasonal_GDP,order = c(8,1,5), seasonal = list(order=c(0,0,1), period=12))
## Warning in arima(deseasonal_GDP, order = c(8, 1, 5), seasonal = list(order
## = c(0, : possible convergence problem: optim gave code = 1
fit_w_seasonality3
## 
## Call:
## arima(x = deseasonal_GDP, order = c(8, 1, 5), seasonal = list(order = c(0, 0, 
##     1), period = 12))
## 
## Coefficients:
##          ar1     ar2      ar3     ar4      ar5     ar6     ar7      ar8
##       0.7845  0.5484  -0.6352  0.0461  -0.6726  0.5098  0.7706  -0.9503
## s.e.  0.0843  0.1010   0.2014  0.1101   0.1721  0.1447  0.1342   0.2729
##           ma1      ma2     ma3     ma4     ma5    sma1
##       -0.7261  -0.9143  0.7954  0.0696  0.0850  0.9391
## s.e.   0.3262   0.2716  0.3524  0.2514  0.1956  1.9185
## 
## sigma^2 estimated as 0.2867:  log likelihood = -53.71,  aic = 137.42
seas_fcast3 <- forecast(fit_w_seasonality3, h=5)
plot(seas_fcast3)

tsdisplay(residuals(fit_w_seasonality3), lag.max = 18, main = "(8,1,5)(0,0,1)[12] Residuos del Modelo Estacional Ajustando el Autoregresor")

#Probando el modelo estacional con Promedio Móvil y Autoregresor ajustado
hold2 <- window(ts(deseasonal_GDP), start=50)
fit2_no_hold2 <- arima(ts(deseasonal_GDP[-c(50:55)]), order=c(8,1,5),seasonal = list(order=c(0,0,1), period=12))
## Warning in arima(ts(deseasonal_GDP[-c(50:55)]), order = c(8, 1, 5),
## seasonal = list(order = c(0, : possible convergence problem: optim gave
## code = 1
fcast2_no_hold2 <- forecast(fit2_no_hold2, h=5)
plot(fcast2_no_hold2)
lines(ts(deseasonal_GDP))