svm-financas
Igor Alexandre Clemente de Morais
May 17, 2017
Support Vector Regression GARCH (SVR-GARCH)
A discussão sobre a caracterização da volatiliade de um ativo via modelos GARCH é vasta e envolve formulações das mais variadas. Cada modelo com uma proposta de investigação de fatos estilizados que ajudam a entender de forma mais precisa as oscilações ao longo do tempo. Porém, o sucesso da maioria desses modelos em caracterizar o passado por vezes não é reproduzido nas previsões. Nesse caso, modelos de redes neurais acabam desempenhando um papel melhor, vide Taylor(2000) e Dunis e Huang(2002), especialmente com a proposta de Vapnik(1995,1997) de usar o algoritmo SVM para estimar a volatilidade, no que é denominado de modelo SVR-GARCH. O que o SVR faz é aproximar uma função que possa representar a relação entre as características (variáveis independentes x) com o resultado (variável dependente y). Quanto mais características tivermos, ou seja, quanto maior for a dimensão do espaço onde estão os “inputs”, mais perfeita será essa aproximação. A conecção entre o espaço das “m” características (x) e o espaço do output (y) é feita por (w), que é um conjunto de pesos lineares \[w=[{w_1},{w_2},...,{w_m}]^T\]. Aqui, a função não-linear \[\Phi{(x)}=[\Phi_1{(x)},....,\Phi_m{(x)}]^T\], com j=1,2,…,m representa as características. Assim, existe uma função de decisão \[f(x)\] que conecta o espaço das características com o espaço dos outputs: \[f(x)=\sum^m_{j=1}w_j\phi_j(x)+b\] Aqui não conhecemos os coeficientes \[w_j\] e também o parâmetro b, denominado de “viés”. Os mesmos devem ser estimados a partir dos dados disponíveis para que seja possível encontrar a função de decisão \[f(x)\]. E podemos usar uma função kernel para determinar \[\phi(x)\]. Dois pontos a destacar aqui. O primeiro é que não existe uma função kernel que seja definida para todos os problemas, podendo a mesma variar de acordo com as características dos dados, ver a discussão sobre machine learning que contém as diversas opções de função kernel a serem utilizadas. O segundo ponto importante nesse processo de estimação de \[f(x)\] é determinar se existe ou não uma função perda, que penaliza erros maiores que determiado valor de \[\epsilon\]. Vapnik(1995) propos uma função perda: \[L_{\epsilon}(x,y,{f(x)})={|y-fx)|}-\epsilon \ \ \ \ for\ \ \ \ |y-f(x)|>=\epsilon \\ =0 \ \ \ \ \ de \ outra \ forma \] onde é possível ver que erros menores que \[\epsilon\] não são penalizados. Mas, também temos como alternativa especificar que \[f(x)\] pode ser encontrado sem uma função perda, o que é o mesmo que determinar que \[\epsilon=0\], ou seja, que o erro é igual a zero. No primeiro exercício abaixo, seguimos a proposta de Pérez-Cruz et al(2003), que usa esse algoritmo para mostrar como pode ser feita a estimativa de um modelo de volatilidade do tipo GARCH(1,1). A ideia aqui é formular a equação da volatilidade diretamente e usar os conceitos de SVR - Support Vector Regression para encontrar os parâmetros de interesse. Duas vantagens podem ser apontadas nesse procedimento. A primeira diz respeito ao fato de não ser necessário impor, a priori, qual seria a distribuição dos resíduos em um modelo de Máxima Verossimilhança. A segunda é em relação a previsão, onde é possível identificar que o algoritmo svm possui melhores resultados na previsão da volatilidade do que os modelos da família GARCH. Vale destacar que o processo de otimização no algoritmo SVM é com base em uma função quadrática e com restrição linear, resultando em um único valor mínimo, sem que se tenha diversas outras soluções locais. Essa propriedade é interessante para problemas no qual temos a presença de não-liearidade nos dados, gerando ganhos sobre estimativas que façam uso de uma função de densidade. ### Estimando um SVR para volatilidade Como é sabido, o modelo GARCH(1,1) pode ser dado a partir de \[y_t=\phi+\sigma_t\epsilon_t \\ \sigma_t=w+\alpha{y^2_t}+\beta{\sigma_{t-1}} \] onde os valores de \[\sigma_t\] podem ser determinados a partir do conhecimento de \[\epsilon_t\], supondo que o mesmo segue uma função de distribuição. A proposta de usar o SVM para estimar \[\sigma_t\] parte apenas do conhecimento da variável em si. Em uma estimativa que usamos a função de verossimilhança os parâemtros são determinados de forma conjunta. A proposta de usar o SVR-GARCH difere do GARCH paramétrico pois não especificamos como seria o processo gerador de dados, ou seja, não precisamos determinar qual é a função de distribuição. Nesse caso, precisamos apenas determinar qual é a variável de resposta e os atributos, também conhecidos em econometria como variáveis independentes. Na proposta de Pérez-Cruz et al(2003) elevamos os retornos ao quadrado, que seriam utilizados posteriormente na equação GARCH(1,1). Já a estimativa da volatilidade pode ser feita de duas formas. Primeiro usamos o desvio padrão histórico dos retornos: \[\sigma_t=\sqrt\frac{\sum_{i=1}^{n}(x_i-\bar{x})^2}{n-1}\] escolhendo uma janela para sua estimativa. A seguir, podemos determinar que a volatilidade é dada pela média ponderada dos retornos ao quadrado: \[\sigma_t=\sqrt\frac{\sum_{i=0}^{n}(y^2_{t-k})}{n}\] A partir de então temos a série dos retornos ao quadrado e seu desvio-padrão, e podemos especificar a fórmula para a volatilidade: \[\sigma_t=w+\alpha{y^2_{t-1}}+\beta{\sigma_{t-1}}\] Aqui é usado o processo conhecido como “feedforward SVR” que possui a mesma estrutura AR, porém, resulta em uma menor capacidade de representar a memória de longo prazo, ver Chen et al.(2009). Para exemplificar o uso do SVR-GARCH(1,1) tal como proposto por Pérez-Cruz et al.(2003), vamos usar os dados de fechamento do índice bovespa de 17-05-2006 a 16-05-2017 em um total de 2.717 dias.
dados<-read.table("svrgarch.csv",header=T,sep=",")Como a serie está em número índice, vamos encontrar os retornos e depois o retorno ao quadrado.
logr<-as.numeric(na.omit(diff(log(dados[,2])))*100)
logr2<-logr^2Agora precisamos da estimativa de volatilidade que, como primeiro exemplo, será dada pelo desvio-padrão histórico dos retornos. Usaremos a função runsd() que está no pacote caTools para fazer esse cálculo. Aqui basta informar o vetor de dados e o tamanho da janela móvel que será utilizada. Outra alternativa de estimar a volatilidade é fazer a média ponderada dos retornos ao quadrado. Nesse caso, vamos usar uma janela de 5 dias usando a função runmean().
require(caTools) #usado para encontrar a janela móvel do desvio-padrão
sdr<-runsd(logr,40)
sdr1<-runmean(logr2,5)par(mfrow=c(1,2))
plot(sdr,main="Desvio-padrão dos retornos",ylab="Desvio-padrão",type="l",col="red")
plot(sdr1,main="Retornos ao quadrado",ylab="Retornos ao quadrado",type="l")
Como mostra a equação do SVR-GARCH(1,1), devemos usar o componente AR(1) do desvio-padrão. Assim, devemos aplicar a primeira defasagem tanto no retorno ao quadrado, quanto nas duas formas propostas para avaliar o desvio-padrão, a estimativa do desvio-padrão e a média dos retornos ao quadrado. Para aplicar o lag, usamos a função abaixo, também disponibilizada em ^(https://heuristically.wordpress.com/2012/10/29/lag-function-for-data-frames/).
lagpad <- function(x, k) {
if (!is.vector(x))
stop('x tem que ser um vector')
if (!is.numeric(x))
stop('x tem que ser numeric')
if (!is.numeric(k))
stop('k tem que ser numeric')
if (1 != length(k))
stop('k tem que ser um número único')
c(rep(NA, k), x)[1 : length(x)]
}logr2_lag<-lagpad(logr2,1)
sdr_lag<-lagpad(sdr,1)
sdr1_lag<-lagpad(sdr1,1)Por fim, vamos juntar todas as séries de interesse em um único banco de dados. Para esse exemplo, precisamos da série do desvio-padrão (sdr), do quadrado dos retornos (logr2_lag) e da primeira defasagem do desvio-padrão (sdr_lag), para a primeira estimativa, e (sdr1,sdr1_lag)
yx1<-na.omit(cbind(sdr,logr2_lag,sdr_lag))
yx2<-na.omit(cbind(sdr1,logr2_lag,sdr1_lag))Para estimar o SVR-GARCH(1,1) usamos a função svm() que está no pacote e1071.
São vários os tipos de kernel que podemos usar bem como também podemos ter diferentes parâmetros de escolha para o processo de otimização. Vamos começar com uma kernel linear e especificando um cost=8. No caso de um modelo do tipo SVR - Support Vector Regression, essa penalidade se aplica ao erro de previsão \[e_i=y_i-w^Tx_i-b\]. Assim, quando esse erro for menor que um determinado “epsilon”, em termos absolutos, não haverá penalização. Por outro lado, para erros maiores que um determinado “epsilon”, é aplicada uma penalização linear.
m1<-svm(sdr~logr2_lag+sdr_lag,yx1,type="eps-regression",kernel="linear",cost=8)
m2<-svm(sdr1~logr2_lag+sdr1_lag,yx2,type="eps-regression",kernel="linear",cost=8)Podemos encontrar os parâmetros do nosso modelo em duas partes. Na primeira determinamos o w, que são os parâmetros das variáveis utilizadas. No nosso exemplo temos duas variáveis, a que está relacionada ao retorno ao quadrado logr2 e a que está relacionada a primeira defasagem da volatilidade sdr1, respectivamente alfa e beta no GARCH. A persistência encontrada fica em 0,99, bastante elevada.
w1<-t(m1$coefs)%*%m1$SV
w1## logr2_lag sdr_lag
## [1,] 0.006896803 0.9927004persistencia1<-w1[,1]+w1[,2]
persistencia1## logr2_lag
## 0.9995972w2<-t(m2$coefs)%*%m2$SV
w2## logr2_lag sdr1_lag
## [1,] 0.02930204 0.902421persistencia2<-w2[,1]+w2[,2]
persistencia2## logr2_lag
## 0.931723Para concluir o modelo vamos determinar a constante.
b1<-m1$rho
b1## [1] 0.001869153b2<-m2$rho
b2## [1] 0.0102084Assim, a nossa equação final da volatilidade para o modelo SVR-GARCH(1,1) com estimativa da volatilidade pelo desvio-padrão, teria o seguinte formato. \[\sigma_t=0,0018-0,006{y^2_t}+0,9927{\sigma_{t-1}}\] Podemos agora usar os parametros para estimar a volatilidade realizada e comparar a mesma com os dados originais. Veja que, pelo gráfico, a estimativa ficou muito próxima do desvio-padrão.
sdr_prev1<-as.data.frame(predict(m1,yx1))
sdr_prev2<-as.data.frame(predict(m2,yx2))plot(sdr,main="volatilidade estimada e desvio-padrão",ylab="Vol",type="l",col="red")
lines(sdr_prev1,type="l",col="blue")
legend(2000,5,c("Vol","desvio"),lty=c(1,1),col=c("blue","red"))
E podemos confirmar esse baixo erro usando a estatística de comparação dos modelos RMSE via função rmse() que está no pacote “hydroGOF”. Veja que temos um erro bem baixo.
require(hydroGOF)
RMSE_linear1<-rmse(sdr_prev1,yx1[,1])
RMSE_linear1## predict(m1, yx1)
## 0.05471977require(hydroGOF)
RMSE_linear2<-rmse(sdr_prev2,yx2[,1])
RMSE_linear2## predict(m2, yx2)
## 2.03761Um dos problemas ao usar o algoritmo SVM é determinar o valor de epsilon e do parâmetro “cost”, que irá penalizar os elevados erros da estimativa. A alternativa é testar diferentes modelos até determinar qual seria o melhor. Porém, há uma forma mais fácil de fazer isso que é usar o comando tune(), que permite selecionar o melhor modelo svm a partir da análise de sensibilidade dos parâmetros epsilon e cost, ou ainda gamma e demais.
m1_tune<-tune(svm,sdr~.,data=as.data.frame(yx1),ranges=list(epsilon=seq(0,1,0.5),cost=1:20))plot(m1_tune,main="Resultados para os modelos svm")
Agora podemos determinar qual seria o melhor modelo estimado. Veja que a escolha é pela kernel “radial” com o parâmetro cost=5 e gamma=0.5.
r1<-m1_tune$best.model
r1##
## Call:
## best.tune(method = svm, train.x = sdr ~ ., data = as.data.frame(yx1),
## ranges = list(epsilon = seq(0, 1, 0.5), cost = 1:20))
##
##
## Parameters:
## SVM-Type: eps-regression
## SVM-Kernel: radial
## cost: 3
## gamma: 0.5
## epsilon: 0
##
##
## Number of Support Vectors: 2715Fazer a previsão da volatilidade.
sdr_prev11<-predict(r1,yx1)E encontrar o RMSE. Note que o melhor modelo tem um rmse menor, com uma kernel=“radial” e diferentes parâmetros. Mas, veja que o ganho com o rmse é pequeno.
RMSE_best1<-rmse(sdr_prev11,yx1[,1])
RMSE_best1## [1] 0.05366859Podemos identificar os resultados das estimativas, com a combinação dos parâmetros, pode ser vista usando o comando summary(m1_tune). Ao ver o gráfico da comparação da volatilidades, desvio-padrão, modelo linear e com kernel radial.
plot(yx1[,1],col="blue",type="l")
lines(sdr_prev1,col="red")
lines(sdr_prev11,col="black",pch=6)
legend(2000,5,c("Desvio","Vol-linear","Vol-radial"),lty=c(1,1,1),col=c("blue","red","black"))
Podemos repetir os passos acima usando, ao invés do desvio-padrão, os retornos ao quadrado.
SVR-GARCH(2,1)
Mais defasagens para a volatilidade, como por exemplo a formulação GARCH(2,1), mostrada abaixo, pode ser feita naturalmente usando o algoritmo SVM, bastando usar mais lags para o desvio-padrão. \[\sigma_t=w+\alpha{y^2_t}+{\beta_1}{\sigma_{t-1}}+{\beta_2}{\sigma_{t-2}}\] Como primeiro passo determinamos a segunda defasagem para o desvio-padrão.
sdr_lag2<-lagpad(sdr,2)A seguir, fazemos a composição do novo banco de dados com essa nova variável.
yxx<-na.omit(cbind(sdr,logr2_lag,sdr_lag2))E, por fim, usamos o algoritmo para encontrar o melhor modelo.
m3<-tune(svm,sdr~.,data=as.data.frame(yxx),ranges=list(epsilon=seq(0,1,0.5),cost=1:20))r2<-m3$best.model
r2##
## Call:
## best.tune(method = svm, train.x = sdr ~ ., data = as.data.frame(yxx),
## ranges = list(epsilon = seq(0, 1, 0.5), cost = 1:20))
##
##
## Parameters:
## SVM-Type: eps-regression
## SVM-Kernel: radial
## cost: 4
## gamma: 0.5
## epsilon: 0
##
##
## Number of Support Vectors: 2714Fazer a previsão da volatilidade do SVR-GARCH(2,1).
sdr1_prev<-predict(r2,yxx)
RMSE_best3<-rmse(sdr1_prev,yxx[,1])
RMSE_best3## [1] 0.075736plot(yxx[,1],col="blue",type="l")
lines(sdr_prev11,col="red")
lines(sdr1_prev,col="black",pch=6)
legend(2000,5,c("Desvio","Vol-linear","Vol-radial"),lty=c(1,1,1),col=c("blue","red","black"))
SVR-GARCH recorrente
Outra forma de estimar o modelo SVR-GARCH é usando a proposta de Chen et al(2009), com a aplicação do algoritmo SVR de forma recurssiva. Primeiro o mesmo é usado para a equação da média dos retornos um AR(1), para obter os resíduos. A seguir, esses resíduos são elevados ao quadrado para obeter uma estimativa da variância. Nesse primeiro passo temos apenas o componente ARCH da equação, dado pelo resíduo ao quadrado. Assim, fazemos uma primeira estimativa da volatilidade e encontramos os resíduos (w), calculamos sua primeira defasagem (componente GARCH do modelo) e refazemos a estimativa via SVR. De forma geral, o modelo SVR-GARCH(1,1) é dado por: \[y_t=f(y_{t-1})+u_t \\ u^2_t=g(u^2_{t-1},w_{t-1})+w_t\] Como pode ser visto, esse novo procedimento considera a presença da defasagem dos resíduos ao quadrado da equação da média AR(1) e também da variância estimada. Esse processo é repetido a cada passo fazendo a atualização do modelo. Em seu artigo, Chen et al(2009) fazem previsões tanto um passo a frente quanto vários passos a frente para a volatilidade. Os resultados apontados atestam a superioridade do modelo SVR-GARCH em previsão. Como primeiro passo vamos abrir os dados.
bov<-read.table("bovespa.csv",header=T,sep=",")
bov<-ts(bov,start = c(2006,5), frequency = 252 )Remove linhas vazias ou com informação específica, como “null”, por exemplo.
bov<-bov[!apply(bov=="",1,all),]Analisando apenas os dados de preços de fechamento ajustados.
serie<-bov[,6]
logr<-as.numeric(na.omit(diff(log(serie)))*100)Fazemos o lag dos retornos e agregamos com os retornos.
logr1<-lagpad(logr,1)
x<-na.omit(cbind(logr,logr1))
head(x,2)## logr logr1
## [1,] -0.1959228 -1.2720611
## [2,] -3.3304978 -0.1959228O primeiro passo é a estimativa de um modelo AR(1) para os retornos para obter os resíduos: \[r_t=\alpha+\beta{r_{t-1}}+u_t\] Fazemos a estimativa do modelo AR(1) considerando o melhor modelo a partir da função tune() e pegamos os residuos dessa equação.
require(e1071)
ar<-tune(svm,logr~.,data=as.data.frame(x),ranges=list(epsilon=seq(0,1,0.5),cost=1:20))
ar1<-ar$best.model
u<-ar1$residualsA seguir, elevamos os resíduos ao quadrado para serem utilizados na estimativa da volatilidade. Aproveitamos para criar a primeira defasagem dos resíduos ao quadrado. Podemos ver a presença de um dos fatos estilizados mais comuns nesse tipo de exercício que é a existência de cluster de volatilidade, seja em períodos de baixa ou de alta.
u2 <- na.omit(cbind(u, lagpad(u,1)))^2 #eleva os resíduos ao quadrado
colnames(u2)<-c("u2","u2_t1")par(mfrow=c(1,2))
plot(u,type="l",main="Residuos do modelo AR(1)",yla="retornos")
plot(u2[,1],type="l",main="Resíduos ao quadrado",ylab="u2")
Agora podemos usar o algoritmo SVM para estimar a volatilidade. Nesse primeiro passo temos uma equação que contém apenas a volatilidade em t e em t-1. Aqui o termo (w) é o resíduo da estimativa da volatilidade. \[\sigma_t=\alpha+\beta\sigma_{t-1}+w_t\]
svm_sd<-svm(u2~.,data=as.data.frame(u2),gamma=5000,kernel="radial",epsilon=0.0001,cost=1000)Podemos fazer a previsão da volatilidade e plotar seus resultados juntamente com os resídulos do modelo AR(1) ao quadrado obtido no passo 1.
sdr3_prev<-predict(svm_sd,u2)plot(u2[,1],col="blue",type="l")
lines(sdr3_prev,col="red")
legend(1500,100,c("desvio","SVR-ARCH(1) recursivo"), lty=c(1,1), col=c("blue","red"))
Abaixo calculamos o rmse para o modelo acima, que mostra uma baixa aderência da nossa estimativa.
require(hydroGOF)
RMSE0<-rmse(sdr3_prev,u2[,1])
RMSE0## [1] 7.056343Mas, o modelo GARCH não está completo. Após a estimativa anterior, encontramos o valor de w (residuo) do modelo e geramos sua primeira defasagem (abaixo uma outra maneira de mostrar como fazer isso sem ser com a formula anterior).
w<-svm_sd$residuals
w1<-c(0, w[-length(w)])
plot(w1,type="l")
Agora temos os elementos para poder gerar uma equação GARCH(1,1) completa. Assim, agregamos o componente \[w_{t-1}\] (update), como variável de informação adicional para a nova estimativa da volatilidade e teremos a equação: \[\sigma_t=\alpha+\beta\sigma_{t-1}+\gamma{w_{t-1}}+w_t\] Vamos criar o banco de dados novo, com a presença de \[w_{t-1}\].
u2w<- cbind(u2[,1:2], w1)
head(u2w,3)## u2 u2_t1 w1
## 2 11.2346386 0.04703187 0.00000
## 3 2.0981039 11.23463860 10.23036
## 4 0.7777945 2.09810391 1.92992svm_sd2<-svm(u2~.,data=as.data.frame(u2w),gamma=4000,kernel="radial",epsilon=0.0001,cost=500)sdr4_prev<-predict(svm_sd2,u2w)plot(u2w[,1],col="blue",type="l")
lines(sdr4_prev,col="red")
legend(1500,100,c("desvio","SVR-GARCH(1,1) recursivo"), lty=c(1,1), col=c("blue","red"))
A seguir, determinamos o rmse do modelo GARCH, que se mostra bem melhor que o anterior que não considerava a presença de \[w_{t-1}\].
RMSE1<-rmse(sdr4_prev,u2w[,1])
RMSE1## [1] 5.764369Prevendo a volatilidade
Um dos problemas em fazer a previsão da volatilidade é que a mesma não é conhecida o que demanda que seja definida a forma de estimar essa volatilidade ex-post. Há várias propostas na literatura da área, como usar os retornos ao quadrado, ou então a diferença entre os retornos e sua média elevado ao quadrado. Como a média é muito próximo de zero, essas duas estatísticas terão resultados com pouca diferença.
Nosso desafio é criar um algoritmo que permita fazer isso de forma recorrente, fazendo a atualização pelo resíduo w a cada passo. Vamos fazer isso em dois passos. Primeiro separamos uma parte dos dados para treino.
n=60 #total de informações treino<-nrow(u2)-n u2_tr <- u2[1:treino,] #treino head(u2_tr,2)
n=60 eps<-0.05 cost<-0.1 gamma<-1/(2*0.05^2) i <- 1 p_count <- 0
while(p_count < 5){ print(i) #imprime o número de loops #Estima o SVM para u^2 svmr<-svm(x=u_tr[,-1],y=u_tr[,1],scale=FALSE,type=“eps-regression”,kernel=“radial”,gamma=svm_gamma,cost=svm_cost,epsilon=svm_eps) #Test autocorrelation of residuals to lag 1 test<-Box.test(svmr\(residuals,lag = 1,type = "Ljung-Box") p_val <-test\)p.value #p-valor p_count <- ifelse(p_val > 0.05, p_count + 1, 0) #Extract residuals for next estimation step w<-svmr$residuals w<-c(0, w[-length(w)]) #lag 1 u_tr <- cbind(u_tr[,1:2], w) i <- i + 1 }
cv<-c(1,2,39,89) uu<-c(0,cv[-length(cv)]) cbind(cv,uu)
Bibliografia
Chen,S.; Härdle, W.K.; Jeong,K.; Forecasting Volatility with Support Vector Machine-Based GARCH model, Journal of Forecasting, 29(4), 406-433 DOI: 10.1002/for.1134, 2010.
Tang,L-B; Tang, L-X; Sheng, H-Y; Forecasting volatility based on wavelet support vector machine, Expert Systems with Applications: An International Journal, V. 36 (2), 2901-2909, 2009
Bezerra, P.C.S. & Albuquerque, P.H.M. - Volatility forecasting via SVR-GARCH with mixture of Gaussian kernels, Computational Management Science (2017) 14: 179. doi:10.1007/s10287-016-0267-0