Modelo de Bornholdt para simulação de mercados financeiros artificiais

Grupo: Leonardo Barcelos, Luana Bianchi e Rubens Borrasca

OBJETIVO

Modelo de Bornholdt

Para estudar os fenômenos físicos que nos rodeiam, usamos modelos matemáticos para entender como sistemas evoluem com o tempo. No entanto, esses modelos muitas vezes podem ser base para estudos além do sistema físico de interesse. Para um sistema de spins, um dos modelos mais simples que leva em conta apenas a interação de cada spin com seus primeiros vizinhos, é o Modelo de Ising. O hamiltoniano que o descreve é defindo pela segunte expressão:

 ${\mathcal {H}}=\sum _{\langle i,j\rangle }J_{ij}S_{i}^{z}S_{j}^{z}$

onde $J_{ij}$ é o termo que define a força de interação entre os spins i e j e $S^{z}$ a componente Z do spin. Além disso $\langle i,j\rangle$ denota que estamos somando somente sobre os primeiros vizinhos. Para um valor de $J_{ij}>0$ , a interação ferromagnética é favorecida, e os spins vizinhos tendem a se alinhar conjuntamente, formando domínios magnéticos. No entanto, para $J_{ij}<0$ a interação antiferromagnética é favorecida, e os spins tende a se "anti-alinhar".

Simulação de Mercados de Ações

Podemos fazer um paralelo interessante entre os spins e a sua interação com seus vizinhos com traders em mercados finaceiros e a suas interações com outros traders. No caso do Modelo de Ising, podemos interpretar cada spin como um trader, ou agente, em um mercado de ações e a direção da componente Z desse spin como a sua estratégia de compra ou venda. Então neste caso mais simples, para $J_{ij}>0$ , os agentes tenderão a adotar a mesma estratégia de mercado, e isso não descreve bem a realidade dos mercados. Para o caso em que $J_{ij}<0$ , os agentes irão adotar estratégias contrárias aos seus vizinhos, e isso também não nos levará a uma representação realista do mercado financeiro. Para que o modelo possa se aproximar da realidade, devemos incrementar o hamiltoniano do sistema.

Um modelo bastante próximo da realidade dos mercados de ações é o proposto por Lux e Marchesi [referências 12 e 13 bornholdt] que classificam os agentes em duas estratégias: os fundamentalistas e os chartistas. Este modelo reproduz muitas das propriedades observadas em mercados reais, como lei de potências distribuídas dos retornos $ret(t)=ln\left({\frac {p(t)}{p(t-1)}}\right)$ dos preços $p(t)$ e uma alta correlação da volatilidade de preços. Além disso, no nível de estratégias, apresenta o fenômeno em que o número de chartistas (ou traders ruidoso)[talvez tirar essa parte] se correlacionam com fases de grande volatilidade, como visto em mercados reais. Este modelo tem grande sucesso, no entanto sua complexidade é alta. Stefan Bornholdt então propôs um modelo de máxima simplicidade, baseado no modelo de Ising, desenhado para simular a dinâmica de expectativas em sistemas de muitos agentes.

No modelo de Bornholdt, há pelo menos duas forças conflitantes vistas em ações econômicas:

O quê os vizinhos fazem: comumente associado à ação de chartistas;
O quê a minoria faz: comportamento associado a fundamentalistas, traders com maior conhecimento sobre valores fundamentais dos mercados de ações.

Neste modelo, essas duas interações conflitantes são combinadas: as interações entre vizinhos são representadas pelo modelo de Ising mais simples; um acoplamento à minoria como um observável global é introduzido por um acoplamento à magnetização global do sistema de spins. Assim, o hamiltoniano incrementado fica com a seguinte expressão:

 ${\mathcal {H}}=\sum _{\langle i,j\rangle }J_{ij}S_{i}S_{j}-\sum _{i}\alpha C_{i}{\frac {1}{N}}\sum _{j}S_{j}$

onde $\alpha$ é o termo de acomplamento de cada spin com a magnetização da rede $(\alpha >0)$ e $C_{i}$ representa a estratégia com relação à magnetização da rede.

Por simplicidade, assume-se que cada spin é atualizado com uma dinâmica de banho térmico de acordo com:

 $S_{i}(t+1)={\begin{cases}+1,&{\text{com }}p=1/[1+exp(-2\beta h_{i}(t))]\\-1,&{\text{com }}1-p\end{cases}}$

Considerando um modelo com $i=1,...,N$ spins, com as orientações $S_{i}(t)=\pm 1$ , a dinâmica dos spins dependerá do campo local $h_{i}(t)$ :

 $h_{i}(t)=\sum _{j}J_{ij}S_{j}-\alpha C_{i}(t){\frac {1}{N}}\sum _{j}S_{j}(t)$

Neste trabalho, o termo $J_{ij}$ é tomado como uma constante $J$ , e ele é responsável pela indução de ordem ferromagnética local. Em particular, o segundo acomplamento permite casos de ordem ferromagnética local e ordem antiferromagnética global.

Alguns conceitos importantes

Retornos ^[1]

Quando se trata de sistemas financeiros, os estudos se concentram mais no retorno dos ativos do que no preço em si, pois a série temporal dos retornos tem propriedades estatísticas mais interessantes que a série dos preços.

Sendo P(t) o preço de um ativo financeiro no instante t, e P(t-1) o preço do ativo no instante (t-1), o retorno linear do ativo é:

 $r(t)={\frac {P(t)-P(t-1)}{P(t-1)}}$

Reescrevendo esta equação, obtemos que:

 $r(t)+1={\frac {P(t)}{P(t-1)}}$

Aplicando a função logarítmica em ambos os lados da equação, e considerando que:

 $ln(x+1)\approx x$

obtêm-se o retorno logarítmico, que é mais indicado quando se têm ativos voláteis, que possuem uma variação muito alta:

 $r(t)=ln\left({\frac {P(t)}{P(t-1)}}\right)$

Considerando que neste estudo serão comparados retornos de diferentes índices, e também os retornos obtidos através das simulações com o modelo de Bornholdt, é importante normalizar os retornos:

 ${\bar {r}}_{n}={\frac {r_{n}-\langle r\rangle }{\sigma _{r}}}$

em que $\sigma _{r}$ é o desvio padrão da serie de retornos e $\langle r\rangle$ a média.

Distribuição dos Retornos ^[2]

Quando se tem um volume considerável de dados é possível obter a distribuição probabilística deles. Para isso pode-se utilizar a estimação de densidade de Kernel (KDE). Ao observar uma pequena janela de tamanho h em torno de um ponto em análise, pode-se dizer que:

 $P(x-h\leq x\leq x+h)\approx {\frac {1}{N}}\sum _{n}^{N}W(u)$

sendo $W(u)$ uma função kernel e $u$ uma variável tal que:

 $u={\frac {x-x_{n}}{2h}}$

Para este estudo utilizou-se um kernel gaussiano:

 $W(u)={\frac {1}{\sqrt {2\pi }}}e^{-{\frac {1}{2}}u^{2}}$

Este método foi aplicado para as séries de retorno para obter a distribuição deles, utilizando $h=0.5\sigma$ .

Volatilidade

Exemplo de uso da volatilidade como série temporal de um retorno arbitrário.

Uma forma de calcular a volatilidade da série temporal de retornos ao longo do tempo é elevar ao quadrado os valores da série. Deste modo pode-se obter uma variável como a que está ilustrada na figura ao lado.

O interessante em estudar volatilidade de retornos financeiros é que essa variável reflete o quão imprevisível é um determinado ativo. Uma ação com alta volatilidade tende a ter um risco maior de investimento, ao passo que ações com baixa volatilidade geralmente retornam riscos menores, pois seu comportamento acaba sendo mais previsível.

Um fato estilizado financeiro é que a volatilidade das séries temporais de retorno apresentam comportamento sazonal por natureza. Há períodos de alta volatilidade, seguidos por períodos com baixa volatilidade, que então são novamente seguidos por alta volatilidade, e assim adiante. E uma forma de mensurar isto é verificando a presença de clusters de autocorrelação na volatilidade de retornos. Isto é, através da análise da autocorrelação da volatilidade, encontrar bolhas que indiquem as fases destes comportamentos.

Para obter a auto correlação o Teorema de Wiener-Khinchin ^[3] ^[4] foi utilizado, de forma que:

 $A(\tau )={\mathcal {F}}\{|S(\omega )|^{2}\}$

onde $S(\omega )$ é a transformada de Fourier do quadrado dos retornos.

Simulações

Variação do tamanho da grade

Um dos objetivos deste estudo é verificar qual o tamanho de grade gera um resultado que melhor simula um índice financeiro. A escolha de índices financeiros para comparação, ao invés de ações ou commodities, é por causa da instabilidade que ativos financeiros separados têm. Imaginando um cenário fictício onde uma empresa A vende sorvetes enquanto a empresa B vende chocolate quente. A tendência é que haja uma sazonalidade nos 2 ativos, de tal forma que, no verão, as ações de A subam enquanto B desce, e vice-versa no inverno. Porém, enquanto os 2 ativos estão se movimentando, o movimento do mercado representado por estas 2 ações permanece aproximadamente constante. Como índices financeiros são, na verdade, médias de um conjunto grande e diverso de ações presentes no mercado, são melhores para a comparação com o modelo.

Foram escolhidos 4 diferentes tamanhos de grade:

16x16: 256 agentes
32x32: 1024 agentes
50x50: 2500 agentes
100x100: 10000 agentes

Cada uma das simulações se deu sob as seguintes condições:

$\alpha =4$
$J=1$
$T=1,5$
Todos os agentes são racionais, ou seja, eles podem trocar de opinião ao longo da simulação, assumindo $C_{i}=-1$ em alguns momentos e $C_{i}=1$ em outros, de acordo com o grupo que o agente está inserido: maioria ou minoria.

Abaixo é possível observar fotografias do sistema em diferentes instantes para cada uma das simulações:

Situação do sistema de 256 agentes (grade 16 x 16) em 9 instantes diferentes.	Situação do sistema de 1024 agentes (grade 32 x 32) em 9 instantes diferentes.	Situação do sistema de 2500 agentes (grade 50 x 50) em 9 instantes diferentes.	Situação do sistema de 10000 agentes (grade 100 x 100) em 9 instantes diferentes.

Como saída da simulação é obtida uma série temporal dos valores de magnetização, que neste modelo significam o preço de um ativo presente no mercado. Utilizando a série de preços (magnetização), foi obtido, então, a série dos retornos logarítmicos normalizados, através do método explicado na seção dos Retornos. Para comparar o modelo com dados reais, obteve-se as series temporais do preço de índices do mercado financeiro: Ibovespa, S&P 500, Dow Jones e NASDAQ. Esses dados foram obtidos com a biblioteca Pandas Datareader do Python, que permite obter dados de páginas da internet como a Yahoo Finance, que possui a série temporal dos preços de vários ativos financeiros. Assim como nas simulações, a partir da série de preços dos índices foram calculadas a série temporal dos retornos logarítmicos normalizados.

Na figura abaixo pode-se observar os retornos de cada simulação e índice financeiro:

Retornos logarítmicos normalizados obtidos para as simulações com diferentes números de agentes (tamanho da grade) na primeira coluna, e na segunda coluna os retornos logarítmicos normalizados para os índices financeiros Ibovespa, S&P 500, Dow Jones e NASDAQ

Alguns pontos interessantes a se observar:

os retornos dos índices S&P 500, Dow Jones e NASDAQ são muito parecidos, isso de deve ao fato de que todos são índices de mercados dos Estados Unidos, enquanto o Ibovespa é um índice de mercado brasileiro;
os retornos das simulações vão variando mais tempo perto do zero conforme o número de agentes cresce;
conforme o número de agentes no sistema aumenta, diminui-se o número de transições agudas de retorno, representadas pelos picos no gráfico. Isso indica um mercado mais estável.
o objeto da simulação não é retornar uma cópia do observado nos dados reais, mas sim um cenário parecido. O intuito é criar um novo mercado, artificial, com características parecidas ao observado nos mercados reais.

Para ter uma comparação melhor a fim de entender que tamanho de grade simula melhor um mercado financeiro, é importante ver a distribuição dos retornos e a auto correlação das volatilidades. Com base no conteúdo da seção Distribuição dos Retornos foi gerado uma curva de densidade de probabilidade para cada série de retorno, que podem ser observadas na figura abaixo:

Distribuição dos retornos das simulações (linha cheia) e dos índices financeiros (linha tracejada). Percebe-se que a distribuição dos retornos das simulações que mais se aproximam da distribuição dos índices financeiros são as que possuem um menor número de agentes.

Desta figura pode-se também observar alguns pontos:

quanto maior o número de agentes, mais longe a distribuição dos retornos das simulações ficam da distribuição dos retornos dos índices de mercado;
as caudas das distribuições dos retornos das simulações vão ficando mais pesadas conforme o número de agentes diminui, o que indica mais uma vez que as simulações com menor número de agentes possuem distribuição dos retornos mais parecidas com a do mercado financeiro, que também possui uma cauda pesada;
para números muito grandes de agentes, há uma distribuição de retorno de cauda mais leve, o que implica em muito menos situações de alto retorno, um mercado mais estável e consequentemente irreal.

Na figura abaixo estão presentes as auto correlações dos retornos quadrados, ou seja, da volatilidade dos retornos:

Auto correlação das volatilidades dos retornos das simulações (primeira coluna) e do retorno dos índices financeiros (segunda coluna). Em ambas as colunas, a escala dos eixos y e x dos gráficos é logarítmica.

Desta imagem são notados alguns aspectos como:

o comportamento da auto correlação das volatilidades das simulações com mais agentes são mais parecidas com as dos índices financeiros;
embora a simulação com 256 agentes tenha a distribuição dos retornos que mais se aproximou, analisando a sua auto correlação da volatilidade se observa que para este número de agentes os retornos são muito voláteis, e portanto não descrevem tão bem um sistema financeiro;
mesmo que as simulações com maior número de agentes apresente a auto correlação das volatilidades mais similar as do mercado financeiro, a partir da simulação com 1024 agentes, nota-se a presença de clusters de volatilidade, comprovando o fato estilizado de que existe memória do retorno quadrado.

Considerando os pontos apresentados, percebe-se que grades de um tamanho grande não descrevem o mercado financeiro da melhor forma, e da mesma forma, mesmo que aparente descrever bem devido a distribuição de retornos semelhantes a do mercado financeiro, a simulação de 256 agentes é muito volátil se comparada com um mercado real. Com isto, considera-se que a grade de tamanho 32x32 melhor representa um mercado de ações, pois possui uma distribuição dos retornos que não difere tanto da distribuição dos retornos dos índices financeiros, e ainda apresenta o fato estilizado da memória da volatilidade.

Variação da opinião dos agentes

O parâmetro $C_{i}$ , como dito anteriormente, indica a opinião de um agente $i$ presente na rede. Esta opinião está relacionada a seguir (ou não) o comportamento da maioria dos outros agentes presentes no sistema.

Por consequência, pode-se resumir as possíveis opiniões que um agente tem na rede em 3 ramos diferentes:

Opinião 1: $C_{i}=+1$ . Agente que opõe sua opinião à da maioria dos outros agentes no sistema, chamado de fundamentalista. Recebem este nome pois se apoiam no princípio fundamental de oferta e demanda da economia, o qual diz que, com demanda maior, a oferta é menor, enquanto que para demandas menores, a oferta é maior. Logo, opinando diferente da maioria, garante-se um maior retorno.
Opinião 2: $C_{i}=-1$ . Agente que iguala sua opinião à da maioria dos outros agentes no sistema, chamado de chartista. Recebe este nome do termo em inglês chart (gráfico). São agentes que sempre analisam, através de gráficos, as ações em alta, para comprá-las, e as em baixa, para vendê-las. Por isso, seguem a maioria do sistema.
Opinião 3: $C_{i}=c\in \{-1;+1\}$ . Agente que não possui estratégia, e a cada passo da simulação joga aleatoriamente do lado dos chartistas ou fundamentalistas, sem raciocínio. É chamado de completamente irracional.

Abaixo, estão gráficos que mostram o comportamento do sistema composto inteiramente por cada tipo de opinião considerada na rede.

Simulação com 100% da rede contendo agentes de Opinião 1	Simulação com 100% da rede contendo agentes de Opinião 2	Simulação com 100% da rede contendo agentes de Opinião 3

Para o sistema composto apenas pela opinião 1, nota-se uma volatilidade muito grande na série temporal. Como os agentes de opinião 1 sempre jogam contra a maioria, quando o sistema atinge um estado de preços grande o suficiente, por exemplo, ocorre uma "debandada" dos agentes, que passam a opinar contra. Desta forma, o preço abaixa rapidamente, até chegar num estado suficientemente baixo, onde os agentes passam a opinar a favor.

Conclusões

Programas

Código para obter a série temporal do preço dos índices financeiros

import pandas_datareader as pdr
from datetime import datetime

ibov = pdr.get_data_yahoo(symbols='^BVSP',start=datetime(1995,1,1),end=datetime(2021,1,1))
SP500 =  pdr.get_data_yahoo(symbols='^GSPC',start=datetime(1995,1,1),end=datetime(2021,1,1))
DJ = pdr.get_data_yahoo(symbols='^DJI',start=datetime(1995,1,1),end=datetime(2021,1,1))
Nasdaq = pdr.get_data_yahoo(symbols='^IXIC',start=datetime(1995,1,1),end=datetime(2021,1,1))

Código para calcular os retornos e para os normalizar

def ret(x):
    N = len(x)

    y = []
    for i in range(N-1):
        r = np.log(x[i+1])-np.log(x[i])
        y.append(r)

    return y

def normalize(x):
    
    N = len(x)    
    y = []
    
    for i in range(N):
        n = x[i] - np.mean(x)
        n = n/np.std(x)
        
        y.append(n)
        
    return np.array(y)

Código que realiza a estimativa de densidade kernel

def gauss(x,mean,std_dev):
    u = (x - mean) / std_dev
    c = 1 / (np.sqrt(2 * np.pi))
    return c * np.exp(- 0.5 * u ** 2)

def kde(x,kernel="gauss",bw=0.1,n_points=1500):
    kernel_options = ["gauss"]
    data = np.array(x)
    x_kde = np.linspace(np.min(data)-bw,np.max(data)+bw,n_points)
    n = data.shape[0]     #Number of rows
    m = x_kde.shape[0]    #Number of columns
    kde_i = []
    if kernel == kernel_options[0]:
        for x in data:
            kde_i.append(gauss(x_kde,x,bw))
    else:
        print("Kernel not found!")
        print("Kernel options are:")
        for k in kernel_options:
            print(" - " + k)
        return np.nan
          
    kde_i = np.array(kde_i).reshape(n,m)
    
    kde = np.array([np.sum(kde_i[:,i]) for i in np.arange(m)])
    kde_norm = kde / np.sum(kde)
    return x_kde,kde_norm

Código para calcular a auto correlação das volatilidades

from scipy.signal import get_window
from scipy.fft import rfft, rfftfreq, irfft

def fft_calculation(y,window="parzen"):
    w = get_window(window=window, Nx=len(y))
    f = y * w
    freqs = np.fft.rfftfreq(len(y))
    return freqs,np.abs(rfft(f))

def autocorr_calculation(y):
    ift = irfft(np.abs(y) ** 2)
    ift_norm = np.abs(ift) / np.abs(ift).max()
    return ift_norm

w_ibov, s_ibov = fft_calculation(np.array(normalize(r_ibov))**2)
c_ibov = autocorr_calculation(s_ibov)

w_sp500, s_sp500 = fft_calculation(np.array(normalize(r_SP500))**2)
c_sp500 = autocorr_calculation(s_sp500)

w_DJ, s_DJ = fft_calculation(np.array(normalize(r_DJ))**2)
c_DJ = autocorr_calculation(s_DJ)

w_nasdaq, s_nasdaq = fft_calculation(np.array(normalize(r_nasdaq))**2)
c_nasdaq = autocorr_calculation(s_nasdaq)

w100, s100 = fft_calculation(np.array(normalize(R100))**2)
c100 = autocorr_calculation(s100)

w50, s50 = fft_calculation(np.array(normalize(R50))**2)
c50 = autocorr_calculation(s50)

w32, s32 = fft_calculation(np.array(normalize(R32))**2)
c32 = autocorr_calculation(s32)

w16, s16 = fft_calculation(np.array(normalize(R16))**2)
c16 = autocorr_calculation(s16)

Código para gerar a figura dos retornos

fig, ax = plt.subplots(4,2,figsize=(10,7))

ax[0][0].plot(np.arange(len(R16)),R16,'lightsteelblue',label='# agentes = 256')
ax[1][0].plot(np.arange(len(R32)),R32,'cornflowerblue',label='# agentes = 1024')
ax[2][0].plot(np.arange(len(R50)),R50,'blue',label='# agentes = 2500')
ax[3][0].plot(np.arange(len(R100)),R100,'midnightblue',label='# agentes = 10000')

ax[0][1].plot(np.arange(len(r_ibov)),normalize(r_ibov),'pink',label='Ibovespa')
ax[1][1].plot(np.arange(len(r_SP500)),normalize(r_SP500),'palevioletred',label='S&P500')
ax[2][1].plot(np.arange(len(r_DJ)),normalize(r_DJ),'mediumvioletred',label='Dow Jones')
ax[3][1].plot(np.arange(len(r_nasdaq)),normalize(r_nasdaq),'purple',label='NASDAQ')

# Setting labels & titles

fig.suptitle('Retornos(t)',fontsize=14)
fig.text(0.5,0, 't', ha='center', fontsize=12)
fig.text(0.28,0.93, 'Bornholdt', ha='center', fontsize=12)
fig.text(0.77,0.93, 'Índices', ha='center', fontsize=12)
fig.text(0, 0.5, 'r(t)', va='center', rotation='vertical',fontsize=12)


for aa in ax: 
    for a in aa:
        a.xaxis.set_major_locator(plt.MaxNLocator(8))
        a.yaxis.set_major_locator(plt.MaxNLocator(5))
        a.legend(loc='upper left')

fig.tight_layout()
plt.show()

Código para gerar a figura da distribuição dos retornos

y100 = np.array(R100)
y50 = np.array(R50)
y32 = np.array(R32)
y16 = np.array(R16)

kde100 = kde(y100,"gauss",np.std(y100)/3,1500)
kde50 = kde(y50,"gauss",np.std(y50)/3,1500)
kde32 = kde(y32,"gauss",np.std(y32)/3,1500)
kde16 = kde(y16,"gauss",np.std(y16)/3,1500)

y_ibov = np.array(normalize(r_ibov))
y_sp500 = np.array(normalize(r_SP500))
y_dj = np.array(normalize(r_DJ))
y_nasdaq = np.array(normalize(r_nasdaq))

kde_ibov = kde(y_ibov,"gauss",np.std(y_ibov)/3,1500)
kde_sp500= kde(y_sp500,"gauss",np.std(y_sp500)/3,1500)
kde_dj = kde(y_dj,"gauss",np.std(y_dj)/3,1500)
kde_nasdaq = kde(y_nasdaq,"gauss",np.std(y_nasdaq)/3,1500)

fig, ax = plt.subplots(figsize=(10,6))


ax.set_xlabel('Retorno normalizado',fontsize=12)
ax.set_ylabel('Densidade de probabilidade',fontsize=12)
fig.suptitle('Distribuição do retorno (normalizado)',fontsize=14)

plt.plot(kde100[0],kde100[1],label='# agentes 10000')
plt.plot(kde50[0],kde50[1],label='# agentes 2500')
plt.plot(kde32[0],kde32[1],label='# agentes 1024')
plt.plot(kde16[0],kde16[1],label='16')

plt.plot(kde_ibov[0],kde_ibov[1],'--',label='Ibovespa')
plt.plot(kde_sp500[0],kde_sp500[1],'--', label='S&P500')
plt.plot(kde_dj[0],kde_dj[1],'--', label='Dow Jones')
plt.plot(kde_nasdaq[0],kde_nasdaq[1],'--', label='NASDAQ')


ax.set_xlim(-3,3)
plt.legend()
fig.tight_layout()

Código para gerar a figura da auto correlação das volatilidades

fig, ax = plt.subplots(4,2,figsize=(10,7))

ax[0][0].plot(c16[:int(len(c16)/2)],'lightsteelblue',label='# agentes = 256')
ax[1][0].plot(c32[:int(len(c32)/2)],'cornflowerblue',label='# agentes = 1024')
ax[2][0].plot(c50[:int(len(c50)/2)],'blue',label='# agentes = 2500')
ax[3][0].plot(c100[:int(len(c100)/2)],'midnightblue',label='# agentes = 10000')

ax[0][1].plot(c_ibov[:int(len(c_ibov)/2)],'pink',label='Ibovespa')
ax[1][1].plot(c_sp500[:int(len(c_sp500)/2)],'palevioletred',label='S&P500')
ax[2][1].plot(c_DJ[:int(len(c_DJ)/2)],'mediumvioletred',label='Dow Jones')
ax[3][1].plot(c_nasdaq[:int(len(c_nasdaq)/2)],'purple',label='NASDAQ')

# Setting labels & titles
#ax[2].set_xlabel('Data',fontsize=12)
#ax[1].set_ylabel('Retorno',fontsize=12)
fig.suptitle('Autocorrelação das volatilidades',fontsize=14)
fig.text(0.5,0, '$\\tau$', ha='center', fontsize=12)
fig.text(0.28,0.93, 'Bornholdt', ha='center', fontsize=12)
fig.text(0.77,0.93, 'Indíces', ha='center', fontsize=12)
fig.text(0, 0.5, '$A(\\tau)$', va='center', rotation='vertical',fontsize=12)



for aa in ax: 
    for a in aa:
        a.xaxis.set_major_locator(plt.MaxNLocator(8))
        a.yaxis.set_major_locator(plt.MaxNLocator(5))
        a.legend(loc='lower left')
        a.set_yscale('log')
        a.set_xscale('log')
        

fig.tight_layout()
plt.show()

Referências

[ret-1] Retornos. Portal Action.

[kde-2] Kernel Density Estimation. Wikipedia

[wiener1-3] Wiener-Khinchin Theorem. Wolfram Math World

[wiener2-4] Wiener–Khinchin theorem. Wikipedia

[1]

[2]

[3]

[4]

Modelo de Bornholdt para simulação de mercados financeiros artificiais