Heterocedasticidade e modelos GARCH

A volatilidade é uma característica de algo que é inconstante, que muda com facilidade. Por exemplo, uma aplicação financeira é volátil quando seu preço sofre alterações frequentes, sendo um investimento mais arriscado. Ela pode ser uma volatilidade histórica (estimada em um período no passado), implícita (gerada por eventos, como eleições presidenciais ou divulgação de balanços) ou realizada (revisão de uma volatilidade prevista).

Uma quantificação da volatilidade pode ser obtida através do cálculo do desvio padrão dos valores de uma série temporal. O desvio padrão é a raiz quadrada da variância, que é o somatório do quadrado das diferenças entre cada valor da amostra em relação a média, dividido pelo número de amostras menos 1. A variância também pode ser calculada com relação a um modelo (em vez da média). Caso o valor da variância não seja igual para todas as observações, é dito que o modelo apresenta heterocedasticidade – “scedasticidade” significa dispersão – caso contrário, apresenta homocedasticidade. No caso de uma série temporal, a variância varia em função do tempo.

No caso de ações, é muito comum séries de retornos ao invés do uso da série dos preços dos ativos. O retorno é definido pela razão entre: a diferença do preço (p) no instante t pelo preço no instante anterior (t-1) e o preço nesse mesmo instante anterior. Além de ser uma variável mais interessante para o investidor, essa série costuma apresentar estacionariedade, fraca dependência linear e não linear, caudas pesadas na distribuição e excesso de curtose, além de comportamento heterocedástico condicional.

\(r_t=log\left ( \frac{p_t}{p_{t-1}} \right )\)

\(p_t=p_{t-1}e^{r_t}\)

A heterocedasticidade pode apresentar aglomerados de volatilidade e choques. O choque é uma mudança de sentido no comportamento da variável, podendo ser positivo ou negativo, e também é chamado de inovação.

Modelos GARCH

Uma distribuição de probabilidade pode ser resumida por características conhecidas como momentos da função de distribuição. A média e a variância são as medidas-resumo (ou momentos) mais frequentemente utilizadas, conhecidos como momentos de primeira e de segunda ordem. Às vezes, é preciso também considerar os momentos de ordem superior. Os momentos de terceira e quarta ordem estão associados ao cálculo de assimetria e de curtose, respectivamente.

Modelos financeiros e macroeconômicos mais antigos tendem a enfatizar somente o primeiro momento condicional (ou seja, considerando uma probabilidade condicional). Neles, as dependências temporais de ordens superiores são tratadas como perturbações aleatórias em seus momentos incondicionais. Essas dependências expressam a existência de aglomerações de volatilidade na série, ou seja, a alternância de períodos de baixa volatilidade com períodos de alta volatilidade. Novos modelos foram desenvolvidos considerando a distinção de uso entre os momentos de segunda ordem condicionais e não-condicionais – a matriz de covariância condicional depende dos estados passados, variando no tempo. Isso é importante para as séries financeiras, nas quais é comum a ocorrência de valores extremos.

A matriz de covariância não-condicional das variáveis de interesse pode ser invariante no tempo, mas no caso da condicional, ela depende de estados passados. Além disso, as séries com valores extremos costumam não ter uma distribuição normal.

O primeiro modelo desenvolvido para estimação da volatilidade foi proposto por Engle (1982) denominado ARCH (Autoregressive Conditional Heterocedasticity). Sua ordem q define a dependência temporal de sua variância – por exemplo, para q = 1, sua variância condicional depende apenas do instante anterior. Na construção de modelos ARCH, um primeiro passo é tentar ajustar um modelo ARIMA, para remover a correlação serial na série, se esta existir.

Em 1986, Bollerslev propôs uma generalização dos modelos ARCH, então conhecidos como GARCH (generalized ARCH). Neles, é adicionada uma componente referente à variância condicional nos instantes anteriores, permitindo descrever a volatilidade com menos parâmetros que um ARCH. Ou seja, a variância do erro está relacionada com os termos de erro elevados ao quadrado em vários períodos anteriores, em vez de somente o período anterior. Um modelo GARCH(p,q) tem sua equação da variância descrita como:

\(\epsilon_t=\sigma_t\varepsilon_t\)

\(\sigma^2_t=\omega+\displaystyle\sum_{i=1}^{p}\alpha_i\epsilon^2_t+\displaystyle\sum_{j=1}^{q}\beta_j\sigma^2_{t-j}\)

onde \(\sigma^2_t\) é a variância no instante t, \(\omega>0\) é a constante e \(\varepsilon\) é o erro i.i.d. Os coeficientes \(\alpha\) e \(\beta\) estão relacionados à persistência de inovação a curto e longo prazo, respectivamente – a persistência da volatilidade indica o quanto uma nova informação pode causar um choque.

Considerando-se a semelhança entre os modelos ARMA(max(p,q),q) e GARCH(p,q), as funções de autocorrelação (FAC) e de autocorrelação parcial (FACP) devem sugerir se a série é heterocedástica. A FAC dá a ordem máxima da auto-regressão do GARCH e a FACP dá a ordem p das “médias móveis” do GARCH. Se p < max(p,q), seus valores são exatamente as ordens dos modelos, caso contrário, q é encontrado por tentativa e erro. Isso vale para modelos simétricos na variância.

Testes de heterocedasticidade

Para verificar se a série apresenta heterocedasticidade condicional, pode-se aplicar o teste de Ljung-Box (uma generalização do teste de Box-Pierce) na série temporal da variável. Nele, a hipótese nula (H0) é de que os resíduos são independentes e identicamentes distribuídos (i.i.d), ou seja, não existe diferença significativa entre a amostra e a população de referência.

Também existem os testes de Engle , que na verdade aplica os testes de Portmanteau Q e de Lagrange. O teste de Portmanteau Q é semelhante ao teste de Ljung-Box, ao examinar se os quadrados de resíduos são uma sequência de ruído branco, enquanto que o teste de multiplicador de Lagrange serve para ajustar um modelo de regressão linear para os resíduos quadrados e examinar se o modelo ajustado é significativo. Neles, a hipótese nula é que os resíduos quadrados são uma sequência de ruído branco, ou seja, os resíduos são homocedásticos.

Aplicando o teste de Ljung-Box para analisar a autocorrelação da variável, quando aceita-se a hipótese de que os retornos são i.i.d., então não é preciso utilizar um modelo ARIMA para ajustar os dados, podendo-se ajustar um modelo (G)ARCH diretamente. Utilizando o mesmo teste no quadrado dos resíduos, o modelo é apropriado para descrever a dependência linear entre os sucessivos retornos no caso de um p-valor alto.

Modelagem e previsão

Geralmente, utiliza-se o critério de Akaike, Bayes ou Hannan-Quinn para escolher o tipo de modelo GARCH. Para sua estimação, é usado o método de máxima verossimilhança. Pode-se considerar a hipótese de distribuição normal, mas no caso de leptocurtose (achatamento da distribuição devido a valores além de três desvios padrão), sugere-se o uso da distribuição t-student ou de erro generalizado (GED). A obtenção dos parâmetros iniciais pode ser feita por Mínimos Quadrados Ordinários (MQO) ou ARIMA.

Em modelos econométricos com variância não-condicional, a incerteza do erro de predição é crescente com o horizonte de previsão. Na presença de modelos GARCH, o intervalo de confiança exige a avaliação das variâncias condicionais do erro futuro. Em modelos com propriedades assintóticas conhecidas, como o GARCH(1,1), já existem funções definidas.

Assim como o ARIMA, a previsão é feita para um passo no futuro, que é incorporado à série para fazer o passo seguinte, e assim sucessivamente. De modo geral, a volatilidade é mais previsível no longo prazo, se comparado com o curto prazo.

GARCH Multivariado

Dentre as formulações mais comuns, está o modelo de correlação condicional dinâmica de Engle e Sheppard (2002), ou DCC. Ele é uma generalização do modelo de correlação condicional constante (CCC), permitindo que os modelos multivariados sejam estimados a vários passos: primeiro, estimar os modelos individualmente assumindo uma estrutura GARCH para então a correlação dinâmica a partir dos resíduos padronizados.

A cada passo no tempo, a covariância assintótica dos parâmetros do estimador da correlação precisa ser ajustada. A matriz de variância-covariância, proposta por Bollerslev & Wooldridge (1992), de cada modelo permanece consistente.

Para implementar a estimação por máxima verossimilhança, primeiro decompõe-se a densidade conjunta da amostra recursivamente, como um produto das densidades condicionais. Assumindo-se normalidade condicional para as inovações do modelo GARCH(p,q), a função de verossimilhança condicional e derivada da normal pode ser escrita.

Implementação em R

O pacote aTSA (Alternative Time Series Analysis) possui a função “arch.test”, que executa os testes Portmanteau Q e Lagrange Multiplier para a hipótese nula de que os resíduos de um modelo ARIMA são homocedásticos. Esse modelo vem da função “estimate”, disponível no mesmo pacote.

O exemplo a seguir abre um arquivo CSV com as linhas no formato “date,value” (YYYY-MM-DD,float) de frequência diária e cria uma série temporal. Essa série é usada para criar o modelo ARIMA (parâmetros estimados previamente), que serve de entrada para o teste:

require(aTSA)

# Read CSV and create time series object
df = read.csv('data.csv', as.is = TRUE)
ano_ini = as.numeric(substr(df$date[1], start = 1, stop = 4))
mes_ini = as.numeric(substr(df$date[1], start = 6, stop = 7))
dia_ini = as.numeric(substr(df$date[1], start = 9, stop = 11))
freq = 365.25
vec_start = c(ano_ini,mes_ini,dia_ini)
serie_model = ts(df[, name_col], freq = freq, start = vec_start)

# Calculate model
model = estimate(vaz_obs, p = 5, d = 1, q = 0, PDQ = spdq)
# Do ARCH test
arch_test = arch.test(model, output = TRUE)
print(arch_test)

Os pacotes “rugarch” e “mgarch” rel=”noopener” target=”_blank”>mgarch” possuem implementações para usar modelos GARCH. Como o DCC é multivariado, é preciso construir um objeto XTS contendo a variável endógena e as exógenas (nesse caso, duas séries temporais no total).

require(xts)
require(rugarch)
require(rmgarch)

# Create XTS array
dates_vec_time = as.POSIXct(df$date, format = "%Y-%m-%d")
endo_obs_xts = xts(x = endo_obs, order.by = dates_vec_time)
exo_obs_xts = xts(x = exo_obs, order.by = dates_vec_time)
data_xts = merge(endo_obs_xts, exo_obs_xts, all = TRUE)

# Create GARCH(1,1)
garch11 = ugarchspec(mean.model = list(armaOrder = c(1,1), include.mean = FALSE), variance.model = list(garchOrder = c(1,1)), distribution.model = "norm")
# DCC specification
dcc.garch11 = dccspec(uspec = multispec(replicate(2,garch11)), dccOrder = c(1,1), distribution = "mvnorm")
# Model
dcc.fit11 = dccfit(dcc.garch11, data = data_xts)
# Forecast
dcc.focast = dccforecast(dcc.fit11, n.ahead = 1, n.roll = 0)
print(dcc.focast)

O parâmetro “n.ahead” define o número de previsões a frente. Quando n.roll > 0 e n.ahead = 1, essa é uma previsão de rolagem pura com base nos dados de saída disponíveis fornecidos na chamada para a rotina de ajuste. Também é assumido que, se algum regressor externo for passado para a rotina de ajuste, eles contêm valores suficientes para cobrir o período out.sample, para que possam ser usados nesse cenário de previsão.

Implementação em python

Para estimar a ordem do modelo GARCH(p,q), primeiro pode ser construído um gráfico de autocorrelação dos resíduos quadrados (variância, que é a grandeza modelada). Ele resume graficamente a força de um relacionamento com uma observação em uma série temporal com observações em etapas anteriores. Os valores que estivem fora do intervalo de confiança (95%, por padrão) definem os potenciais valores dos parâmetros, pois é muito provável que cada um seja uma correlação e não um acaso estatístico. A autocorrelação parcial no atraso k é a correlação resultante após a remoção do efeito de qualquer correlação devido aos termos em intervalos mais curtos.

def acf(self, df):
	"""
	Plot ACF and PACF for estimate p and q
	"""
	import numpy as np
	from matplotlib import pyplot
	from statsmodels.graphics.tsaplots import plot_acf
	from statsmodels.graphics.tsaplots import plot_pacf
	
	# Plot values
	#pyplot.plot(df['y'])
	#pyplot.show()
	# Calculate retorno
	retorno = np.diff(np.log(df['y'].values))
	#pyplot.plot(retorno)
	#pyplot.show()
	
	# Plot ACF
	plot_acf(retorno)
	pyplot.show()
	
	# Plot PACF
	plot_pacf(retorno, lags=50)
	pyplot.show()

O pacote “arch” deve ser instalado no python pelo canal “bashtage”, através do comando “conda install -c bashtage arch”. Ele possui alguns dados de valores de ações que podem ser usados.

Observando-se a descrição do modelo gerado, a hipótese nula é de que os coeficientes (gerados pela matriz das variâncias) são significativos e existe heterocedasticidade. Assim, se o p-valor for próximo de zero, é verificada essa hipótese e existe uma interferência de inovação.

É possível gerar previsões para processos GARCH (p, q) padrão usando qualquer um dos três métodos de geração de previsão: analítico, baseado em simulação e baseado em Bootstrap.

O script a seguir calcula um modelo GARCH(1,1). Primeiro, calcula o retorno sobre uma série com preços das ações das 500 maiores empresas da bolsa de Nova York. Depois, é definido e calculado o modelo para a série para então calcular a previsão em um horizonte de 5 valores.

def garch(self, df):
	from arch import arch_model
	import arch.data.sp500
	
	data = arch.data.sp500.load()
	market = data['Adj Close']
	returns = 100 * market.pct_change().dropna()

	# Alternative returns calculation
	#returns = np.diff(np.log(market.values))
	
	am = arch_model(returns, vol = Garch', p = 1, o = 0, q = 1, dist = 'Normal')
	res = am.fit(disp = 'off')
	print(res)
	forecasts = res.forecast(horizon = 5)
	# Alternative forecast using simulation
	#forecasts = res.forecast(horizon = 15, method = 'simulation', simulations = 1000, rng = None)
	
	print(forecasts.mean.iloc[-3:])
	print(forecasts.residual_variance.iloc[-3:])
	print(forecasts.variance.iloc[-3:])

As previsões começam com a especificação do modelo e estimativa de parâmetros. Elas estão contidas em um objeto ARCHModelForecast que possui 4 atributos: mean, residual_variance, variance e simulations (nula para modelo usando “analytics”).

Cada coluna “h.n” da saída da previsão corresponde a n passos a frente. A saída é alinhada para que a coluna “Date” represente os dados finais usados para gerar a previsão, de modo que “h.1” na linha “2018-12-31” seja a previsão de um passo à frente feita usando dados até 31 de dezembro de 2018 inclusive, por exemplo.

Fontes

• Econometria Básica, Damodar N. Gujarati e Dawn C. Porter, McGraw Hill (2011)
• Econometria de Séries Temporais, Rodrigo de Losso da Silveira Bueno, Fundação Getúlio Vargas (2008)
• Portal Action – Modelos GARCH
• Uma análise empírica da volatilidade do retorno de commodities agrícolas utilizando modelos ARCH: os casos do café e da soja, Washington Santos da Silva; Thelma Sáfadi; Luiz Gonzaga de Castro Júnior. Rev. Econ. Sociol. Rural vol.43 no.1 Brasília Jan/Mar. 2005
• Modelagem de Séries Temporais com Python – SciPy-SP (slides de palestra com exemplos de códigos)
• How to Model Volatility with ARCH and GARCH for Time Series Forecasting in Python (tutorial de modelagem)
• Documentação arch – Forecasting (exemplos de modelagem e previsão)