Oilprice Forecast using GMLSS R Package
Leandro Correa
11/05/2023
Problema de negócio
Para ilustrar o uso do pacote gamlss, vamos abordar um problema de negócio interessante: a previsão probabilística do preço do petróleo com antecedência de um dia.
Nesse cenário, a empresa solicitou que fornecêssemos não apenas uma previsão pontual, mas sim uma previsão probabilística, ou seja, uma distribuição de probabilidade para os preços do petróleo.
O pacote gamlss é uma ferramenta poderosa para modelagem estatística e previsão, e pode ser encontrado em (https://cran.r-project.org/web/packages/gamlss/gamlss.pdf).
O conjunto de dados que temos disponível é o ‘oil dataset’ dentro do pacote gamlss, onde a variável ‘OILPRICE’ representa nossa variável de resposta.
Nosso objetivo é desenvolver um modelo preditivo para a variável ‘OILPRICE’ e fornecer a previsão para um dia adiante. Durante o tutorial, iremos explicar como selecionar o modelo adequado, as métricas de diagnóstico utilizadas para avaliar o desempenho do modelo, entre outros aspectos relevantes.
Carregando os pacotes necessários para análise
packages <- c("tidyverse","ggrepel","caret","knitr","kableExtra",
"splines","reshape2","PerformanceAnalytics","correlation","see",
"ggraph","nortest","rgl","car","ggside","olsrr",
"jtools","ggstance","GGally","cowplot","Rcpp",
"gamlss","gamlss.add","gamlss.dist")
sapply(packages, require, character.only = TRUE)## tidyverse ggrepel caret
## TRUE TRUE TRUE
## knitr kableExtra splines
## TRUE TRUE TRUE
## reshape2 PerformanceAnalytics correlation
## TRUE TRUE TRUE
## see ggraph nortest
## TRUE TRUE TRUE
## rgl car ggside
## TRUE TRUE TRUE
## olsrr jtools ggstance
## TRUE TRUE TRUE
## GGally cowplot Rcpp
## TRUE TRUE TRUE
## gamlss gamlss.add gamlss.dist
## TRUE TRUE TRUEFunções e variáveis globais
# It's a metric used to evaluate the accuracy of prediction models or estimates.
mape <- function(y, yhat) {
mean(abs((y - yhat) / y)) * 100
}
remove_corr <- function(df, pred_variables, cut_off = 0.85) {
# Calculate correlation matrix
cor_matrix <- cor(df %>% dplyr::select(-one_of(pred_variables)))
# Find indices of variables with correlation > cut_off
high_cor <- findCorrelation(cor_matrix, cutoff = cut_off)
# Remove highly correlated columns
df_filtered <- df[,-high_cor]
df_filtered <- cbind(df_filtered, df[,pred_variables])
return(df_filtered)
}
# Number of observations that will be part of the validation set.
# 25 days is the approximate number of days in a month, excluding weekends.
VALID_DATASET_SIZE <- 25
# K is a required penalty for the GAIC function.
K_SET <- c(2, 2.5, 3, 3.5, 4)Carregamento dos dados
O conjunto de dados “oilprice” contém informações sobre o valor do preço do petróleo em diferentes períodos. Abaixo a descrição de cada variável que compõe o dataset.
Para obter mais informações sobre o conjunto de dados “oilprice”, consulte Stasinopoulos et al., 2015, p. 413.
# Loading oil price Data
data(oil)
Pré-processamento
Como premissa para o desenvolvimento do problema, consideramos o cenário em que conseguimos obter todas as informações que desejamos utilizar como variáveis explicativas com um atraso de um dia (d-1).
Por esse motivo, criamos um segundo conjunto de dados (df) considerando todas as variáveis explicativas com um deslocamento de período de tempo d-1. Dessa forma, cada valor de OILPRICE é representado pelos valores das variáveis independentes no período de tempo anterior. Esse mesmo princípio também é aplicado à variável respLAG, que já representa a variável OILPRICE com um dia de atraso. Após este procedimento temos a variável de resposta OILPRICE contendo os valores atuais do preço do petroleo, OILPRICE.1 passa a representar o valor de OILPRICE a partir do dia d-1, LagResp a partir do dia d-2 e todas as outras variáveis, conforme descrito na seção ‘carregamento dos dados’, também apresentam um atraso igual a d-1.
# One-time unit shift process
df <- data.frame(OILPRICE = oil[-1, "OILPRICE"], oil[-nrow(oil),])
# Checking data consistency
glimpse(df) ## Rows: 999
## Columns: 26
## $ OILPRICE <dbl> 4.633077, 4.634049, 4.646312, 4.631520, 4.627616, 4.635214,…
## $ OILPRICE.1 <dbl> 4.640923, 4.633077, 4.634049, 4.646312, 4.631520, 4.627616,…
## $ CL2_log <dbl> 4.636475, 4.645352, 4.637831, 4.638315, 4.650526, 4.635893,…
## $ CL3_log <dbl> 4.641116, 4.649857, 4.642466, 4.642562, 4.654722, 4.640344,…
## $ CL4_log <dbl> 4.644968, 4.653484, 4.646312, 4.646120, 4.658047, 4.644199,…
## $ CL5_log <dbl> 4.648038, 4.656338, 4.649665, 4.648708, 4.660321, 4.646984,…
## $ CL6_log <dbl> 4.649761, 4.657858, 4.651672, 4.649952, 4.661267, 4.648421,…
## $ CL7_log <dbl> 4.650908, 4.658711, 4.653103, 4.650621, 4.661740, 4.649378,…
## $ CL8_log <dbl> 4.651863, 4.659564, 4.654341, 4.651099, 4.662117, 4.650335,…
## $ CL9_log <dbl> 4.652340, 4.660037, 4.655293, 4.651577, 4.662401, 4.651099,…
## $ CL10_log <dbl> 4.651672, 4.659374, 4.655007, 4.651004, 4.661645, 4.650813,…
## $ CL11_log <dbl> 4.650621, 4.657952, 4.653865, 4.649570, 4.659942, 4.649570,…
## $ CL12_log <dbl> 4.648613, 4.655578, 4.651672, 4.647080, 4.656908, 4.646984,…
## $ CL13_log <dbl> 4.646120, 4.652531, 4.648804, 4.644102, 4.653484, 4.643910,…
## $ CL14_log <dbl> 4.643236, 4.649187, 4.645736, 4.640923, 4.649761, 4.640537,…
## $ CL15_log <dbl> 4.639765, 4.645256, 4.642081, 4.636960, 4.645544, 4.636475,…
## $ BDIY_log <dbl> 6.850126, 6.850126, 6.879356, 6.882437, 6.896694, 6.913737,…
## $ SPX_log <dbl> 7.221624, 7.235309, 7.222756, 7.222252, 7.237620, 7.233556,…
## $ DX1_log <dbl> 4.386554, 4.379762, 4.387449, 4.383675, 4.382364, 4.382951,…
## $ GC1_log <dbl> 7.413367, 7.419680, 7.418481, 7.410347, 7.399704, 7.404888,…
## $ HO1_log <dbl> 1.136197, 1.152564, 1.155182, 1.136743, 1.139946, 1.137256,…
## $ USCI_log <dbl> 4.108412, 4.120986, 4.115127, 4.103965, 4.107096, 4.097672,…
## $ GNR_log <dbl> 3.917806, 3.942552, 3.923952, 3.925531, 3.941970, 3.936325,…
## $ SHCOMP_log <dbl> 7.744539, 7.762536, 7.766061, 7.765158, 7.755763, 7.775213,…
## $ FTSE_log <dbl> 8.636699, 8.650062, 8.639729, 8.642292, 8.659907, 8.656137,…
## $ respLAG <dbl> 4.631812, 4.640923, 4.633077, 4.634049, 4.646312, 4.631520,…Identificando e removendo features com alta correlação
Antes de realizar qualquer modelagem utilizando GLM (Generalized Linear Models), é recomendado identificar e remover variáveis explicativas com alta correlação. Abaixo listamos alguns dos principais motivos:
Redundância de informações: Variáveis altamente correlacionadas podem fornecer informações semelhantes ao modelo, o que pode levar a uma redundância desnecessária. Isso pode aumentar a complexidade do modelo sem adicionar informações significativas.
Instabilidade dos coeficientes: Quando duas variáveis estão altamente correlacionadas, é possível que o modelo atribua pesos excessivos a uma delas, enquanto a outra pode ter um peso próximo a zero. Isso pode levar a instabilidade nos coeficientes estimados e dificultar a interpretação dos resultados.
Violação de pressupostos: A alta correlação entre variáveis pode violar pressupostos de independência dos erros, o que pode levar a resultados incorretos ou enviesados na modelagem.
A função chart.Correlation mostra a distribuição do dataset df, bem como o índice de correlação entre as variáveis explicativas do problema.
chart.Correlation((df[2:ncol(df)]), histogram = TRUE)
Como é possível observar no gráfico acima, há várias variáveis que apresentam um alto grau de correlação (basta verificar que, por exemplo, quanto maior o tamanho do número no triângulo superior da matriz, maior é o grau de correlação entre as variáveis). Com o objetivo de mitigar esse problema, realizamos a análise da correlação de Pearson entre as variáveis explicativas e aplicamos um limiar de 85% para remover as variáveis que apresentavam uma correlação alta.
# Selecting variables that do not have high correlation
df <- remove_corr(df, pred_variables = c("OILPRICE","OILPRICE.1"))
exp_variables <- colnames(df)[colnames(df) %nin% "OILPRICE"]
summary(df) ## CL14_log CL15_log HO1_log USCI_log
## Min. :3.594 Min. :3.603 Min. :-0.1442 Min. :3.650
## 1st Qu.:4.111 1st Qu.:4.117 1st Qu.: 0.6232 1st Qu.:3.838
## Median :4.497 Median :4.494 Median : 1.0547 Median :4.021
## Mean :4.332 Mean :4.331 Mean : 0.8607 Mean :3.962
## 3rd Qu.:4.527 3rd Qu.:4.522 3rd Qu.: 1.1013 3rd Qu.:4.070
## Max. :4.654 Max. :4.649 Max. : 1.1877 Max. :4.148
## GNR_log FTSE_log respLAG OILPRICE
## Min. :3.317 Min. :8.568 Min. :3.266 Min. :3.266
## 1st Qu.:3.787 1st Qu.:8.715 1st Qu.:3.967 1st Qu.:3.966
## Median :3.869 Median :8.778 Median :4.517 Median :4.517
## Mean :3.818 Mean :8.760 Mean :4.311 Mean :4.309
## 3rd Qu.:3.920 3rd Qu.:8.813 3rd Qu.:4.580 3rd Qu.:4.580
## Max. :3.985 Max. :8.868 Max. :4.705 Max. :4.705
## OILPRICE.1
## Min. :3.266
## 1st Qu.:3.966
## Median :4.517
## Mean :4.310
## 3rd Qu.:4.580
## Max. :4.705Ao final, as variáveis: CL14_log, CL15_log, HO1_log, USCI_log, GNR_log, FTSE_log, respLAG e OILPRICE.1 apresentaram um índice de correlação dentro do threshold estipulado e foram selecionadas para a fase de construção do modelo.
Particionando os dados entre treino e teste
Na fase de particionamento dos dados, selecionamos as últimas 25 linhas do conjunto de dados df (ou seja, os últimos 25 dias) para o nosso conjunto de validação. Além disso, selecionamos também a última linha do conjunto de dados oil para nossa previsão probabilistica do dia d+1.
As 25 observações selecionadas foram utilizadas somente para simular o comportamento do modelo considerando uma previsão dentro de um intervalo de aproximadamente um mês (excluindo os finais de semana) e não para aferir qualidade de ajuste do modelo aos dados.
É importante ressaltar que, para modelos determinísticos, como no caso de modelos GLM que são obtidos por meio da otimização de uma função de verossimilhança de densidade de probabilidade, não é necessário dividir os dados em conjuntos de treinamento e teste. Ao fazer uma estratificação treino e teste com um percentual muito alto em modelos GLM, como no caso do padrão 70-30 comumente aplicado, pode ocorrer viés nos parâmetros identificados como betas da função.
# Set the size of training and validation sets
train_size <- nrow(df) - VALID_DATASET_SIZE
# Create the training and validation sets
train_data <- df[1:train_size, ]
valid_data <- df[(train_size+ 1):nrow(df), ]
# Identifying the explanatory variables of the day ahead forecast
predict_test <- oil[nrow(oil), ] %>%
dplyr::rename(OILPRICE.1 = OILPRICE) %>%
dplyr::select(one_of(exp_variables))Estimando a distribuição de melhor ajuste
Após a fase de pré-processamento dos dados e a identificação das correlações entre as variáveis explicativas, damos início ao o uso do pacote “gamlss”. O primeiro passo nessa etapa é selecionar a família de distribuição adequada para a variável de resposta “OILPRICE”. A função “fitDist()”, disponibilizada pelo pacote, utiliza um conjunto de distribuições predefinidas para ajuste aos dados e escolhe a melhor distribuição com base no critério de informação generalizado de Akaike (GAIC), com uma penalidade padrão de κ = 2. A ordem dos modelos ajustados é exibida de forma crescente, ou seja, do “melhor” para o “pior” (Stasinopoulos et al., 2015, p. 155-156).
Considerando o nosso conjunto de treinamento, a distribuição SHASH (sinh-arcsinh) foi a melhor opção para o ajuste da distribuição à variável OILPRICE.
fitted_dist <- fitDist(OILPRICE, data=train_data, type="realline")fitted_dist$fits## SHASH SHASHo SHASHo2 SEP4 SEP1 SEP3
## -298.421838 -101.038152 -101.038152 -72.287416 -22.708476 -16.189540
## EGB2 JSUo JSU SEP2 ST5 ST2
## -6.849126 2.406725 2.406725 14.111737 50.958931 70.592579
## ST1 ST3 SST SN1 SN2 ST4
## 74.316070 77.007409 80.517703 136.193156 141.139384 276.486015
## GU GT PE2 PE NO TF
## 435.867533 541.256756 662.631593 662.631593 773.554309 775.554309
## TF2 exGAUS LO NET RG
## 775.554309 775.637124 795.872518 841.683004 1129.242605histDist(train_data$OILPRICE,"SHASH",nbins=30, n.cyc=100)
##
## Family: c("SHASH", "Sinh-Arcsinh")
## Fitting method: "nlminb"
##
## Call: gamlssML(formula = train_data$OILPRICE, family = "SHASH")
##
## Mu Coefficients:
## [1] 4.525
## Sigma Coefficients:
## [1] 10.33
## Nu Coefficients:
## [1] 10.63
## Tau Coefficients:
## [1] 12.7
##
## Degrees of Freedom for the fit: 4 Residual Deg. of Freedom 970
## Global Deviance: -306.422
## AIC: -298.422
## SBC: -278.896Procedimento Stepwise
O procedimento stepwise é uma abordagem comumente utilizada na análise de regressão para selecionar um conjunto ótimo de variáveis explicativas (preditoras) para um modelo de regressão. Ele é usado para determinar quais variáveis têm um impacto significativo na variável de resposta e devem ser incluídas no modelo. O pacote gamlss propõe duas estratégias A e B, explicadas com mais detahes no livro. Aplicamos as duas estratégias ao problema para identificar qual das duas gera o melhor ajuste, considerando o índice GAIC como critério de validação.
stepGAICAll.A
A função stepGAICAll.A(), do pacote gamlss, é uma estratégia para selecionar termos aditivos utilizando o GAIC para todos os parâmetros de distribuição, pressupondo uma distribuição de resposta específica (Stasinopoulos et al., 2015, p. 397). Aplicamos a estratégia stepGAICAll.A utilizando a família de distribuição SHASH previamente identificada.
# Set number of cpu's, General Formulation, and k set for GAIC estimation
nC<-detectCores()
FORM<-as.formula(paste("~",paste(paste(paste("(",
names(df)[-1], sep=""),")",sep=""),
collapse="+")))# Creating a null model (without explanatory variables) that will serve as a baseline for the stepGAICAll.A method
m0 <- gamlss(OILPRICE~1,family=SHASH,data=train_data,gd.tol=Inf)# Executing the stepGAICAll.A method considering the explanatory variables of the problem and the SHASH distribution family.
mStraA.SHASH.rs <- stepGAICAll.A(m0, scope=list(lower=~1,
upper=FORM),
ncpus=nC,k=sqrt(log(dim(train_data)[1])))- Observando os parâmetros identificados pelo procedimento stepGAICAll.A.
summary(mStraA.SHASH.rs)## ******************************************************************
## Family: c("SHASH", "Sinh-Arcsinh")
##
## Call: gamlss(formula = OILPRICE ~ OILPRICE.1 + respLAG +
## HO1_log + GNR_log, sigma.formula = ~OILPRICE.1,
## nu.formula = ~USCI_log, tau.formula = ~1, family = SHASH,
## data = train_data, gd.tol = Inf, trace = FALSE)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.089543 0.009735 9.198 < 2e-16 ***
## OILPRICE.1 0.891217 0.032318 27.576 < 2e-16 ***
## respLAG 0.099627 0.032967 3.022 0.00258 **
## HO1_log 0.024257 0.004051 5.988 2.99e-09 ***
## GNR_log -0.018955 0.006954 -2.726 0.00653 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9843 0.4648 4.269 2.16e-05 ***
## OILPRICE.1 -1.5197 0.1050 -14.477 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: log
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.3948 1.0651 3.187 0.001482 **
## USCI_log -0.9402 0.2667 -3.526 0.000442 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2932 0.0558 -5.254 1.83e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## No. of observations in the fit: 974
## Degrees of Freedom for the fit: 10
## Residual Deg. of Freedom: 964
## at cycle: 20
##
## Global Deviance: -5192.908
## AIC: -5172.908
## SBC: -5124.094
## ******************************************************************É importante verificar se os valores de p (p-values) de todos os coeficientes estimados (Mu, Sigma e Tau) estão dentro do intervalo de confiança estatística (< 0.05). É válido lembrar também que o intercepto sempre permanece como parâmetro da função e que os índices AIC e GAIC são parâmetros utilizados para avaliar o ajuste do modelo aos dados da variável de resposta.
Buscando aprimorar o ajuste do modelo, o pacote gamlss oferece algumas funções de suavização (smoothers) com o objetivo de capturar padrões gerais e tendências nos dados, reduzindo o ruído e a variabilidade aleatória. A seguir, reajustamos o mesmo modelo identificado pela função stepGAICAll_A(), incorporando ajustes de parâmetros por meio da função de suavização p-spline pb() (Stasinopoulos et al., 2015, p. 16).
mStraA.SHASH.2.rs = gamlss(formula = OILPRICE ~ pb(OILPRICE.1) + pb(respLAG) +
pb(HO1_log) + pb(GNR_log),
sigma.formula = ~pb(OILPRICE.1),
nu.formula = ~pb(USCI_log),
tau.formula = ~1,
family = SHASH,
data = train_data,
gd.tol = Inf,
trace = FALSE,
method=RS(20))Warning in RS(20): Algorithm RS has not yet convergedPercebe-se após a execução da última célula que existe uma mensagem de alerta indicando que o algoritmo RS não convergiu. Os algoritmos especificados como parâmetros na variável “method” são os estimadores operacionais utilizados para identificar os coeficientes com base na máxima verossimilhança da função de distribuição de probabilidade do parâmetro “family”. O pacote “gamlss” conta com três abordagens diferentes de maximização: RS, CG e MX.
Abaixo, executamos as três abordagens e ao final identificamos o melhor ajuste com base no menor índice GAIC.
- Utilizando o algorithmo CG (veja Stasinopoulos et al., 2015, p. 70)
mStraA.SHASH.2.cg <- gamlss(formula = OILPRICE ~ pb(OILPRICE.1) + pb(respLAG) +
pb(HO1_log) + pb(GNR_log),
sigma.formula = ~pb(OILPRICE.1),
nu.formula = ~pb(USCI_log),
tau.formula = ~1,
family = SHASH,
data = train_data,
gd.tol = Inf,
trace = FALSE,
method=CG(20))## Warning in CG(20): Algorithm CG has not yet converged- Utilizando o algorithmo MX (veja Stasinopoulos et al., 2015, p. 197).
# using the mixed algorithm
mStraA.SHASH.2.mx <- gamlss(formula = OILPRICE ~ pb(OILPRICE.1) + pb(respLAG) +
pb(HO1_log) + pb(GNR_log),
sigma.formula = ~pb(OILPRICE.1),
nu.formula = ~pb(USCI_log),
tau.formula = ~1,
family = SHASH,
data = train_data,
gd.tol = Inf,
trace = FALSE,
method=mixed(10,10))## Warning in CG(n.cyc = n2): Algorithm CG has not yet convergedstepGAICAll.B
Da mesma maneira que aplicamos a estratégia stepGAICAll.A, vamos aplicar a estratégia stepGAICAll.B. Nesta estratégia, todos os parâmetros de distribuição são obrigados a terem os mesmos termos. Ou seja, se um termo X for selecionado, ele será incluído no preditor de todos os parâmetros (Stasinopoulos et al., 2015, p. 399-400). Ao final, pretendemos comparar as estratégias A e B para obter o modelo de melhor ajuste.
m0 <- gamlss(OILPRICE~1,family=SHASH,data=train_data,gd.tol=Inf)mStratB.SHASH.rs <-stepGAICAll.B(m0, scope=list(lower=~1,
upper=FORM),
ncpus=nC,k=sqrt(log(dim(train_data)[1])))Observing the parameters identified by the stepGAICAll.B procedure.
summary(mStratB.SHASH.rs)## ******************************************************************
## Family: c("SHASH", "Sinh-Arcsinh")
##
## Call: gamlss(formula = OILPRICE ~ OILPRICE.1 + respLAG +
## USCI_log, sigma.formula = ~OILPRICE.1 + respLAG +
## USCI_log, nu.formula = ~OILPRICE.1 + respLAG +
## USCI_log, tau.formula = ~OILPRICE.1 + respLAG +
## USCI_log, family = SHASH, data = train_data, gd.tol = Inf,
## trace = FALSE)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.034188 0.008031 -4.257 2.28e-05 ***
## OILPRICE.1 0.925648 0.038238 24.208 < 2e-16 ***
## respLAG 0.075100 0.039012 1.925 0.0545 .
## USCI_log 0.007517 0.005741 1.309 0.1907
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.387 6.075 1.216 0.224
## OILPRICE.1 3.206 7.075 0.453 0.651
## respLAG -3.955 6.838 -0.578 0.563
## USCI_log -2.189 2.533 -0.864 0.388
##
## ------------------------------------------------------------------
## Nu link function: log
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.1987 3.7586 2.181 0.0294 *
## OILPRICE.1 0.8858 4.0552 0.218 0.8271
## respLAG -0.1466 3.9272 -0.037 0.9702
## USCI_log -2.9475 1.5462 -1.906 0.0569 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.569 3.326 0.772 0.440
## OILPRICE.1 3.006 4.194 0.717 0.474
## respLAG -2.659 4.152 -0.640 0.522
## USCI_log -1.092 1.380 -0.792 0.429
##
## ------------------------------------------------------------------
## No. of observations in the fit: 974
## Degrees of Freedom for the fit: 16
## Residual Deg. of Freedom: 958
## at cycle: 20
##
## Global Deviance: -5194.958
## AIC: -5162.958
## SBC: -5084.856
## ******************************************************************Vale ressaltar que alguns dos coeficientes estimados no procedimento stepwise B ficaram fora do intervalo de confiança (p-value < 0.05). No entanto, a remoção desses coeficientes que estão além do nível de significância estatística não teve impacto no ajuste do modelo. Levando em consideração que todos os modelos serão comparados ao final da análise e que faremos uma etapa de diagnóstico para avaliar o ajuste dos dados, decidimos manter a configuração indicada pelo procedimento stepwise B.
- Na célula abaixo Treinamos um novo modelo incluindo um suavizador para os parâmetros selecionados pelo procedimento stepGAICAll.B.
mStratB.SHASH.2.rs = gamlss(formula = OILPRICE ~ pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
sigma.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
nu.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
tau.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
family = SHASH,
data = train_data,
gd.tol = Inf,
trace = FALSE,
method=RS(20))Warning in RS(20): Algorithm RS has not yet converged- Executamos um novo modelo incluindo um suavizador para os parâmetros selecionados pelo procedimento stepGAICAll.B, além de testar os algorithmos CG, MX e RS.
# using the CG algorithm
mStratB.SHASH.2.gc = gamlss(formula = OILPRICE ~ pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
sigma.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
nu.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
tau.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
family = SHASH,
data = train_data,
gd.tol = Inf,
trace = FALSE,
method=CG(20))
# using the mixed algorithm
mStratB.SHASH.2.mx = gamlss(formula = OILPRICE ~ pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
sigma.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
nu.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
tau.formula = ~pb(OILPRICE.1) + pb(respLAG) + pb(USCI_log),
family = SHASH,
data = train_data,
gd.tol = Inf,
trace = FALSE,
method=mixed(10,10))## Warning in CG(n.cyc = n2): Algorithm CG has not yet convergedSelecionando o modelo com o índice GAIC mais baixo
Após a fase de stepwise considerando suavização dos estimadores identificados e diferentes estratégias de máxima verossimilhança, partimos para fase de seleção de modelo com o melhor ajuste a partir do índice GAIC.
O índice GAIC (Generalized Akaike Information Criterion) envolve a combinação do valor do logaritmo da função de verossimilhança do modelo ajustado e um termo de penalização k, que leva em consideração o número de parâmetros estimados no modelo.
A fórmula geral para o cálculo do GAIC é a seguinte:
GAIC = -2 * log(L) + k * p
Onde:
L é o valor da função de verossimilhança do modelo ajustado; p é o número de parâmetros estimados no modelo; e k é o termo de penalização.
Abaixo comparamos todos os modelos estimados a partir do procedimento stepwise (A e B), e analisamos os índices GAIC de cada um desses modelos considerando diferentes termos de penalidade.
# Check, using different k , which model is the best
df_gaic = data_frame()
for(k in K_SET){
df_aux <- GAIC(mStraA.SHASH.rs,
mStraA.SHASH.2.rs,
mStraA.SHASH.2.cg,
mStraA.SHASH.2.mx,
mStratB.SHASH.2.rs,
mStratB.SHASH.2.gc,
mStratB.SHASH.2.mx,
k=k)
rownames(df_aux) <- paste0(rownames(df_aux), "_k", k)
df_gaic <- rbind(df_gaic, df_aux)
}df_gaic %>% arrange(AIC) df AIC
mStratB.SHASH.2.mx_k2 22.72227 -5.185255e+03
mStraA.SHASH.2.mx_k2 14.92935 -5.178015e+03
mStratB.SHASH.2.rs_k2 22.45113 -5.177836e+03
mStraA.SHASH.2.rs_k2 14.95126 -5.177358e+03
mStratB.SHASH.2.mx_k2.5 22.72227 -5.173894e+03
mStraA.SHASH.rs_k2 10.00000 -5.172908e+03
mStraA.SHASH.2.mx_k2.5 14.92935 -5.170551e+03
mStraA.SHASH.2.rs_k2.5 14.95126 -5.169883e+03
mStraA.SHASH.rs_k2.5 10.00000 -5.167908e+03
mStratB.SHASH.2.rs_k2.5 22.45113 -5.166610e+03
mStraA.SHASH.2.mx_k3 14.92935 -5.163086e+03
mStraA.SHASH.rs_k3 10.00000 -5.162908e+03
mStratB.SHASH.2.mx_k3 22.72227 -5.162533e+03
mStraA.SHASH.2.rs_k3 14.95126 -5.162407e+03
mStraA.SHASH.rs_k3.5 10.00000 -5.157908e+03
mStraA.SHASH.2.mx_k3.5 14.92935 -5.155621e+03
mStratB.SHASH.2.rs_k3 22.45113 -5.155385e+03
mStraA.SHASH.2.rs_k3.5 14.95126 -5.154932e+03
mStraA.SHASH.rs_k4 10.00000 -5.152908e+03
mStratB.SHASH.2.mx_k3.5 22.72227 -5.151171e+03
mStraA.SHASH.2.mx_k4 14.92935 -5.148157e+03
mStraA.SHASH.2.rs_k4 14.95126 -5.147456e+03
mStratB.SHASH.2.rs_k3.5 22.45113 -5.144159e+03
mStratB.SHASH.2.mx_k4 22.72227 -5.139810e+03
mStratB.SHASH.2.rs_k4 22.45113 -5.132934e+03
mStraA.SHASH.2.cg_k2 73.75207 -4.507318e+03
mStraA.SHASH.2.cg_k2.5 73.75207 -4.470442e+03
mStraA.SHASH.2.cg_k3 73.75207 -4.433566e+03
mStraA.SHASH.2.cg_k3.5 73.75207 -4.396690e+03
mStraA.SHASH.2.cg_k4 73.75207 -4.359814e+03
mStratB.SHASH.2.gc_k2 111.35570 2.676889e+31
mStratB.SHASH.2.gc_k2.5 111.35570 2.676889e+31
mStratB.SHASH.2.gc_k3 111.35570 2.676889e+31
mStratB.SHASH.2.gc_k3.5 111.35570 2.676889e+31
mStratB.SHASH.2.gc_k4 111.35570 2.676889e+31Como resultado final, considerando os valores de GAIC para os cenários de penalização (2, 2.5, 3, 3.5, 4), o modelo mStratB.SHASH.2.mx foi selecionado como aquele que melhor se ajusta aos dados da variável OILPRICE.
Diagnostic plots
Na fase de diagnóstico, utilizaremos algumas funções gráficas do pacote gamlss para analisar o ajuste dos dados do modelo desenvolvido à variável de resposta OILPRICE. O pacote fornece uma variedade de visualizações dos resultados do ajuste, selecionamos três destas visualizaćões para conduzir esta etapa do pipeline de análises.
Residual plot()
A função plot() produz quatro gráficos para verificar os resíduos quantílicos normalizados (aleatorizados) de um modelo gamlss. A aleatorização é realizada para variáveis de resposta discretas e mistas, bem como para dados de intervalo ou censurados. Os quatro gráficos gerados pela função são:
• resíduos em relação aos valores ajustados do parâmetro µ;
• resíduos em relação a um índice ou covariável especificada;
• uma estimativa de densidade do núcleo dos resíduos;
• QQ-normal plot dos resíduos (Stasinopoulos et al., 2015, p. 422-423).
# Residual plots from the fitted normal model mStratB.SHASH.2.mx
plot(mStratB.SHASH.2.mx)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.007582036
## variance = 1.002913
## coef. of skewness = 0.01901131
## coef. of kurtosis = 3.029456
## Filliben correlation coefficient = 0.9995259
## ******************************************************************O gráfico resultante do nosso modelo considerando a variável de resposta OILPRICE mostrado acima é muito semelhante à Figura 12.3 apresentada em (Stasinopoulos et al., 2015, p. 424). Observe que os resíduos se comportam bem, pois os dois primeiros gráficos: (1) dos resíduos em relação aos valores ajustados de µ e (2) em relação aos índices das observações, mostram uma dispersão aleatória ao redor da linha horizontal em 0, enquanto a estimativa de densidade do núcleo dos resíduos é aproximadamente normal (gráfico 3) e o gráfico de Q-Q normal (4) é aproximadamente linear (com interceptação em 0 e inclinação em 1).
Ao analisar também as estatísticas resumidas dos resíduos quantílicos (output abaixo da figura), vemos que sua média é aproximadamente zero, que sua variância é aproximadamente um, que seu coeficiente de assimetria (baseado em momentos) está próximo de zero, e que seu coeficiente de curtose (baseado em momentos) está próximo de 3. As estatísticas sugerem que os resíduos estão aproximadamente distribuídos de forma normal (r ∼ N(0, 1)), como deveriam ser para um modelo adequado. Além disso, o coeficiente de correlação de Filliben (ou o coeficiente de correlação do gráfico de probabilidade normal) está próximo de 1.
Worm plot()
O gráfico worm plot mostra o quão distantes os resíduos estão de seus valores esperados. As duas curvas elípticas no gráfico representam os intervalos de confiança aproximados de 95%. Se mais de 5% dos pontos estiverem fora das curvas, ou se houver uma clara divergência da linha horizontal, o modelo é inadequado para explicar a variável de resposta. A curva ajustada aos pontos do worm reflete diferentes inadequações no modelo. Por exemplo, se o nível de pontos do gráfico estiver acima de uma linha horizontal na origem, isso indica que a média dos resíduos é muito alta. O capítulo 12 do livro de Stasinopoulos et al., (2015) explica em detalhes como interpretar um worm plot.
# Residual plots from the fitted normal model mStratB.SHASH.2.mx
wp(mStratB.SHASH.2.mx)
Como mostra o gráfico acima, a maioria das observações está dentro da região de “aceitação” entre as duas curvas elípticas, indicando que o modelo se ajusta bem no geral. No entanto, há uma distorção na parte direita do gráfico, onde os pontos apresentam uma maior distância da linha horizontal. Isso sugere que a média dos resíduos no final das observações tende a ser maior, indicando uma perda de ajuste nessa região.
Worm dtop()
A função dtop() verifica visualmente o ajuste de um modelo construindo um intervalo de confiança não paramétrico para uma função de distribuição verdadeira, com base na função de distribuição empírica da amostra. O gráfico de Owen transformado e sem tendência (DTOP) é aplicado aos resíduos quantílicos normalizados ajustados do modelo para verificar sua adequação (Stasinopoulos et al., 2015, p. 433-434).
dtop(mStratB.SHASH.2.mx)
Uma vez que a linha horizontal de cada gráfico DTOP está dentro das faixas de confiança de 95%, concluímos que os resíduos normalizados poderiam ter vindo de uma distribuição normal e, consequentemente, a distribuição assumida da variável de resposta (SHASH) é razoável.
Ajuste do modelo considerando a variável respLAG
Abaixo está uma representação do ajuste do modelo considerando a variável “respLAG”, apenas para visualização do comportamento da função identificada.
plot(OILPRICE~respLAG, col="lightblue", data=train_data)
lines(fitted(mStratB.SHASH.2.mx)[order(train_data$OILPRICE)]~
train_data$OILPRICE[order(train_data$OILPRICE)])
Avaliação
Avaliamos o modelo considerando primeiramente a previsão no intervalo de um mês útil (25 dias), com um nível de significância de 95%.
# Prediction
test_df = valid_data[-1]
test_df$pred <- predict(mStratB.SHASH.2.mx, newdata = test_df, interval = "prediction", level = 0.95)
test_df$OILPRICE <- valid_data$OILPRICE- Calculando o erro percentual médio absoluto das previsões.
# Deviation calculation
test_df['MAPE'] <- apply(test_df, 1, function(x) mape(x['OILPRICE'], x['pred']))
test_df[,c('OILPRICE','pred','MAPE')] %>% arrange(MAPE)## OILPRICE pred MAPE
## 1 3.489819 3.488321 0.04292064
## 2 3.633367 3.635396 0.05586150
## 3 3.673512 3.669716 0.10333868
## 4 3.675288 3.680480 0.14127778
## 5 3.645972 3.640341 0.15444114
## 6 3.542986 3.537017 0.16847352
## 7 3.724488 3.717360 0.19138939
## 8 3.646494 3.639476 0.19244357
## 9 3.674781 3.685187 0.28319515
## 10 3.545586 3.528718 0.47575606
## 11 3.592919 3.611491 0.51690449
## 12 3.644928 3.667383 0.61607125
## 13 3.650658 3.627493 0.63454067
## 14 3.615771 3.643171 0.75778170
## 15 3.597312 3.624739 0.76240926
## 16 3.538057 3.508752 0.82826269
## 17 3.605226 3.639854 0.96049150
## 18 3.683616 3.719066 0.96238599
## 19 3.518980 3.481801 1.05653639
## 20 3.581294 3.535188 1.28741892
## 21 3.645189 3.592827 1.43645414
## 22 3.693867 3.639611 1.46880188
## 23 3.726175 3.670223 1.50158643
## 24 3.649619 3.587658 1.69773554
## 25 3.634951 3.571365 1.74928476- Identificando a margem de desvio da previsão considerando uma probabilidade de 95%.
deviation_margin <- quantile(test_df[,'MAPE'], probs = 0.95)
print(paste0("95% of the predicted values are in a +-",round(deviation_margin,2),"% deviagion margin to real values."))## [1] "95% of the predicted values are in a +-1.66% deviagion margin to real values."Considerando os parâmetros de entrada obtidos no dia d, forneceremos uma previsão probabilística para um dia à frente (d+1), com uma margem de erro de +-1,66%, com 95% de confiança.
Resultados
Resultado da previsão para o próximo dia:
Já identificamos a margem de erro com intervalo de confianća de 95%, agora vamos fazer a previsão do dia d+1. Abaixo os valores das variáveis explicativas que serviram de suporte para previsão.
# input parameters (last line of oil dataset)
print(predict_test)## CL14_log CL15_log HO1_log USCI_log GNR_log FTSE_log respLAG OILPRICE.1
## 1000 3.797061 3.801538 0.169574 3.699325 3.554491 8.728248 3.646494 3.605226- Fazendo a previsão pontual do valor de OILPRICE para d+1.
pred.day_ahead <- predict(mStratB.SHASH.2.mx, newdata= predict_test, type = "response")print(paste0("The predicted value for the d+1 OILPRICE variable is: ", round(pred.day_ahead,6)))## [1] "The predicted value for the d+1 OILPRICE variable is: 3.600918"E por fim vamos gerar uma uma estimativa de densidade de probabilidade com base nos valores preditos para variável OILPRICE considerando o conjunto prediction_test. Isso vai nos dar uma distribuição de probabilidade para os preços do petróleo para a previsão do próximo dia.
pred.prob <-predictAll(mStratB.SHASH.2.mx,newdata=predict_test)A curva de probabilidade é obtida adicionando novos casos ao conjunto de dados e o modelo é reajustado para todos os casos (incluindo obtidos no treinamento e os novos adicionados). É importante observar que a extrapolação (adição de novos casos) deve ser evitada ou tratada com cautela, como discutido na Seção 5.4.1 de Stasinopoulos et al., (2015).
pdf.plot(mu=pred.prob$mu,
sigma=pred.prob$sigma,
nu=pred.prob$nu,
tau=pred.prob$tau,
family=SHASH,
min=3.45,
max= 3.8, step=0.001)
O eixo x do gráfico representa os valores preditos para a variável de resposta (no caso, “y”), enquanto o eixo y representa a função de densidade de probabilidade correspondente a esses valores preditos. A função de densidade de probabilidade fornece informações sobre a probabilidade de ocorrer cada valor de y, levando em consideração a distribuição especificada pelo modelo. Em casos em que a curva da distuibuição é mais pontuda, como no nosso exemplo, isso é indicativo de uma menor dispersão e uma maior concentração dos valores preditos em torno de um valor médio ou esperado, o que indica uma boa estimativa.
Considerações finais
Se você chegou até aqui, parabéns! O objetivo deste post não é apenas fornecer uma visão geral do pacote gamlss do R, mas também comentar e destacar alguns pontos importantes sobre o processo de análise de regressão linear que são frequentemente negligenciados em análises de mercado. Espero que o objetivo tenha sido alcançado e, se tiver algum comentário, correção, ajuste ou sujestão, por favor entre em contato pelo hscleandro@gmail.com.
Bibliografia
STASINOPOULOS, Mikis et al. Flexible regression and smoothing: The GAMLSS packages in R. GAMLSS for Statistical Modelling. GAMLSS for Statistical Modeling, 2015.