Exibir código-fonte para 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:
 <math>\mathcal{H} = \sum_{\langle i, j \rangle} J_{ij} S_i^z S_j^z</math>
onde <math>J_{ij}</math> é o termo que define a força de interação entre os spins '''i''' e '''j''' e <math>S^z</math> a componente Z do spin. Além disso <math>\langle i, j \rangle</math> denota que estamos somando somente sobre os primeiros vizinhos. Para um valor de <math>J_{ij} > 0</math>, a interação ferromagnética é favorecida, e os spins vizinhos tendem a se alinhar conjuntamente, formando domínios magnéticos. No entanto, para <math>J_{ij} < 0</math> a interação antiferromagnética é favorecida, e os spins tende a se "anti-alinhar".

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 financeiro 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 <math>J_{ij} > 0</math>, os agentes tenderão a adotar a mesma estratégia de mercado, e isso não descreve bem a realidade dos mercados.

O Modelo de Bornholdt, se incluid {coisas do artigo}. Assim, o hamiltoniano do sistema fica complementado:
 <math>\mathcal{H} = \sum_{\langle i, j \rangle} J_{ij} S_i S_j - \sum_{i} \alpha C_i(t) \frac{1}{N} \sum_{j} S_j(t)</math>
onde <math>\alpha</math> é o termo de acomplamento do spin com a rede toda <math>(\alpha > 0)</math>

==Alguns conceitos importantes==

===Retornos <ref name='ret'>[http://www.portalaction.com.br/series-temporais/51-retornos Retornos. Portal Action.]</ref>===

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 é:

 <math>r(t) = \frac{P(t) - P(t-1)}{P(t-1)}</math>

Reescrevendo esta equação, obtemos que:

 <math>r(t) + 1= \frac{P(t)}{P(t-1)}</math>

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

 <math>ln(x + 1) \approx x</math>

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

 <math>r(t) = ln\left(\frac{P(t)}{P(t-1)}\right)</math>

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:

 <math>\bar{r}_n = \frac{r_n  - \langle r \rangle}{\sigma_r}</math>

em que <math>\sigma_r</math> é o desvio padrão da serie de retornos e <math>\langle r \rangle</math> a média.

===Distribuição dos Retornos <ref name='kde'>[https://en.wikipedia.org/wiki/Kernel_density_estimation Kernel Density Estimation. Wikipedia]</ref>===
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:

 <math>P(x-h \leq x \leq x+h ) \approx \frac{1}{N}\sum_n^N W(u) </math>

sendo <math>W(u)</math> uma função ''kernel'' e <math>u</math> uma variável tal que:

 <math>u = \frac{x - x_n}{2h}</math>

Para este estudo utilizou-se um ''kernel'' gaussiano:

 <math>W(u) = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}u^2}</math>

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

===Volatilidade===

[[Arquivo:Exemplo volatilidade.png|500px|thumb|right|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 <ref name='wiener1'>[https://mathworld.wolfram.com/Wiener-KhinchinTheorem.html Wiener-Khinchin Theorem. Wolfram Math World]</ref> <ref name='wiener2'>[https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem Wiener–Khinchin theorem. Wikipedia]</ref> foi utilizado, de forma que:

 <math>A(\tau) =\mathcal{F}\{|S(\omega)|^2\}</math>

onde <math>S(\omega)</math> é 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:

*<math>\alpha = 4</math>
*<math>J = 1</math>
*<math>T = 1,5 </math>
*Todos os agentes são racionais, ou seja, eles podem trocar de opinião ao longo da simulação, assumindo <math>C_i = -1</math> em alguns momentos e <math>C_i = 1</math> 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:

{| style="margin-left: auto; margin-right: auto; border: none;"
|-
! [[Arquivo:N16.png|250px|thumb|center|Situação do sistema de 256 agentes (grade 16 x 16) em 9 instantes diferentes.     ]] 
! [[Arquivo:N32.png|250px|thumb|center|Situação do sistema de 1024 agentes (grade 32 x 32) em 9 instantes diferentes.    ]] 
! [[Arquivo:N50.png|250px|thumb|center|Situação do sistema de 2500 agentes (grade 50 x 50) em 9 instantes diferentes.    ]]
! [[Arquivo:N100.jpg|250px|thumb|center|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 [[Modelo_de_Bornholdt_para_simulação_de_mercados_financeiros_artificiais#Retornos|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 [https://pandas-datareader.readthedocs.io/en/latest/ ''Pandas Datareader''] do Python, que permite obter dados de páginas da internet como a [https://finance.yahoo.com/ ''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:

[[Arquivo:Retornos.png|600px|thumb|center|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 [[Modelo_de_Bornholdt_para_simulação_de_mercados_financeiros_artificiais#Distribui.C3.A7.C3.A3o_dos_Retornos_.5B2.5D|distribuição dos retornos]] e a auto correlação das [[Modelo_de_Bornholdt_para_simulação_de_mercados_financeiros_artificiais#Volatilidade|volatilidades]]. Com base no conteúdo da seção [[Modelo_de_Bornholdt_para_simulação_de_mercados_financeiros_artificiais#Distribui.C3.A7.C3.A3o_dos_Retornos_.5B2.5D|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:

[[Arquivo:Dist.png|600px|thumb|center|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;
*

Na figura abaixo estão presentes as auto correlações dos retornos quadrados, ou seja, da [[Modelo_de_Bornholdt_para_simulação_de_mercados_financeiros_artificiais#Volatilidade|volatilidade]] dos retornos:

[[Arquivo:Volatilidade.png|600px|thumb|center|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 auto correlação da volatilidade se observa que para este número de agentes os retornos são muito voláteis;
*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 simulação com 1024 agentes apresenta ''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, não descrevendo muito bem o mercado. Com isto, considera-se que a grade de tamanho 32x32 melhor descreve o mercado, 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 <math>C_i</math>, como dito anteriormente, indica a opinião de um agente <math>i</math> 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: <math>C_i = +1</math>. 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: <math>C_i = -1</math>. 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: <math>C_i = c \in \{-1;+1\}</math>. 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.

{| style="margin-left: auto; margin-right: auto; border: none;"
|-
! [[Arquivo:Opinion1.png|350px|thumb|center|Simulação com 100% da rede contendo agentes de Opinião 1]] 
! [[Arquivo:Opinion2.png|350px|thumb|center|Simulação com 100% da rede contendo agentes de Opinião 2]] 
! [[Arquivo:Opinion3.png|350px|thumb|center|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===
<source lang="python">
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))
</source>

===Código para calcular os retornos e para os normalizar===
<source lang="python">
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)
</source>

===Código para calcular a auto correlação das volatilidades===
<source lang="python">
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)
</source>

===Código para gerar a figura dos retornos===
<source lang="python">
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()
</source>

==Referências==

<references/>