Um modelo de optimização não linear restringida na estimação de um modelo GARCH

Jos´e F´elix Mavungo Um modelo de optimiza¸c˜ao n˜ao linear restringida na estima¸c˜ao de um modelo GARCH Departamento de Matem´atica Faculdade de Ciˆencias da Universidade do Porto 2014 Jos´e F´elix Mavungo Um modelo de optimiza¸c˜ao n˜ao linear restringida na estima¸c˜ao de um modelo GARCH Tese submetida `a Faculdade de Ciˆencias da Universidade do Porto para obten¸c˜ao do grau de Mestre em Engenharia Matem´atica Departamento de Matem´atica Faculdade de Ciˆencias da Universidade do Porto 2014 Para os meus pais, Jos´e Mavungo e Aldenice Pambo. 3 Agradecimentos Tese realizada no ˆambito do mestrado em Engenharia Matem´atica Departamento de Matem´atica Faculdade de Ciˆencias da Universidade do Porto http://www.fc.up.pt/dmat/engmat Agrade¸co a Deus por me ter guardado e protegido durante esses dois longos anos distante da minha gente e terra querida. Agrade¸co a Maria Eduarda da Rocha Pinto Augusto da Silva e Maria do Carmo Vaz de Miranda Guedes, orientadoras desta tese de mestrado, pelo apoio, disponibilidade e partilha de conhecimentos. Agrade¸co ao professor Kevin Sheppard do Departamento de Economia da universidade de Oxford pela ajuda e pelos coment´arios na implementa¸c˜ao do modelo no Matlab. Agrade¸co a organiza¸c˜ao Continuum.org na pessoa de Matt Harward por terem facultado os pacotes Accelerate e IOPro para elabora¸c˜ao deste trabalho. Agrade¸co a Seth Deland da MathWorks pelo seu suporte e valioso tempo despendido na optimiza¸c˜ao das routinas Matlab. Agrade¸co a Sven Leyffer da Argonne National Laboratory pelos seus coment´arios acerca do modelo usado. Agrade¸co aos professores Jos´e Lelo da Universidade 11 de Novembro de Cabinda, No´e Africano Pedro da Universidade Lus´ıadas de Cabinda e Vicente Cuancua da Universidade T´ecnica de Angola (UTANGA) pelos seus valiosos incentivos e me terem mostrado o caminho certo a seguir. Agrade¸co a minha fam´ılia e amigos pelo apoio moral e em especial agrade¸co profundamente aos meus pais que por meio de um sacrif´ıcio e amor incondicional permitiram que eu chegasse at´e aqui. 4 Resumo Os modelos GARCH foram introduzidos nos anos 80 para representar adequadamente caracter´ısticas das s´eries financeiras. Quando se consideram v´arias s´eries financeiras, tem-se um contexto multivariado e o modelo multivariado GARCH, ou mais geral VECH, apresenta problemas de parcim´onia. Foram ent˜ao impostas restri¸c˜oes que levaram `a defini¸c˜ao de outros modelos, como o DVECH e BEKK. Neste trabalho, considera-se uma especifica¸c˜ao do modelo GARCH VECH que ´e mais geral do que o DVECH. A formula¸c˜ao deste modelo ´e feita do ponto de vista de programa¸c˜ao n˜ao linear com restri¸c˜oes, e ´e designada por GARCH NLP. O modelo GARCH NLP ´e usado para modelar os ´ındices SP&500 e FTSE. Verificouse que o modelo GARCH NLP capta melhor a volatilidade dos dados do que os modelos tradicionais. 5 Abstract The GARCH models were introduced in the 80s to adequately represent characteristics of financial series. When we combine financial variables, we get a multivariate context and more general GARCH model, VECH presents problems of parsimony. Were then imposed restrictions that led to the development of other models such as the BEKK and DVECH. This work considers a specification of GARCH VECH, which is more general than the DVECH. The formulation of this model is made from the viewpoint of constrained non linear programming, and is designated by NLP GARCH. The NLP GARCH model is used to model the SP&500 and FTSE indices. It was found that the model NLP GARCH better captures volatility than traditional database models. 6 Conte´udo Resumo 5 Abstract 6 ´ Indice de Tabelas 9 ´ Indice de Figuras 10 1 Introdu¸c˜ao 11 1.1 Alguns fundamentos te´oricos sobre s´eries temporais . . . . . . . . . . . . . . 12 1.2 S´eriestemporais.................................. 12 1.2.1 Estacionariedade ............................. 12 1.2.2 Modelos ARMA(p,q) . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3 Fatos estilizados sobre dados financeiros. . . . . . . . . . . . . . . . . . . . . 15 1.4 Retornosdeativos ................................ 16 1.5 Carteira de retornos e pagamento de dividendos . . . . . . . . . . . . . . . . 17 2 Modelos heteroced´asticos 18 2.1 Modelo GARCH univariado . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.1 Estima¸c˜ao dos parˆametros . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1.2 Testando os efeitos ARCH . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1.3 Outros modelos GARCH . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 Modelos GARCH multivariados . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2.1 ModeloVECH .............................. 23 2.2.2 ModeloDVECH ............................. 24 2.2.3 ModeloBEEK............................... 24 2.2.4 Estacionariedade do modelo VECH . . . . . . . . . . . . . . . . . . . 24 2.2.5 Estima¸c˜ao dos parˆametros . . . . . . . . . . . . . . . . . . . . . . . . 25 3 Aplica¸c˜ao dos modelos GARCH 26 3.1 Propriedades estat´ısticas dos retornos de ´ındices . . . . . . . . . . . . . . . . 26 3.2 Modela¸c˜ao da volatilidade dos retornos de´ındices atrav´es dos modelos GARCH: casounivariado .................................. 30 3.3 Modela¸c˜ao da volatilidade dos retornos de´ındices atrav´es dos modelos GARCH: casomultivariado................................. 32 3.3.1 Avalia¸c˜ao do modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 7 4 Conclus˜oes 40 Bibliografia 41 A Resultados de simula¸c˜ao 44 Referˆencias 44 B Resultados dos modelos DVECH e BEKK. 46 C C´odigos Matlab 50 D C´odigos S-Plus 66 8 Lista de Tabelas 3.1 Medidas de tendˆencia central, de dispers˜ao, de assimetria e de achatamento dosretornos. ................................... 26 3.2 Teste dos efeitos ARCH sobre os retornos. . . . . . . . . . . . . . . . . . . . 30 3.3 Parˆametros estimados, estat´ıstica t-student, erros padr˜ao, fun¸c˜ao de logm´axima-verosimilhan¸ca do modelo GARCH univariado dos retornos dos´ındices S&P500eFTSE.................................. 31 3.4 Resultado do caso bivariado sobre retornos de ´ındices S&P500 e FTSE (valores entre parˆenteses representam os erros padr˜ao), (*) coeficientes significativos com erros padr˜ao robustos e (**) coeficientes significativos com erros padr˜ao robustos e obtidos por aproxima¸c˜ao num´erica. . . . . . . . . . . . . . 34 3.5 Resultado do teste sobre a n˜ao correla¸c˜ao dos res´ıduos padronizados de S&P500 eFTSE. ...................................... 38 3.6 Teste dos efeitos ARCH sobre os retornos de padronizados. . . . . . . . . . . 38 3.7 Teste sobre a normalidade dos retornos padronizados. . . . . . . . . . . . . . 39 A.1 Estat´ıstica descritiva de 10.000 simula¸c˜oes com observa¸c˜oes de 10000 erros N(0,1) ....................................... 44 B.1 Teste de Jarque-Bera sobre a normalidade dos retornos padronizados. . . . . 46 B.2 Teste de Ljung-Box sobre a existˆencia de autocorrela¸c˜ao nos retornos padronizados do modelo DVECH. . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 B.3 Teste de multiplicador de Lagrange (LM) sobre existˆencia de efeitos ARCH nos retornos padronizados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 B.4 Teste de Jarque-Bera sobre a normalidade dos retornos padronizados do modeloBEKK. .................................. 47 B.5 Teste de Ljung-Box sobre a existˆencia de autocorrela¸c˜ao nos retornos padronizadosdomodeloBEKK............................. 48 B.6 Teste de multiplicador de Lagrange sobre existˆencia de efeitos ARCH nos retornos padronizados do modelo BEKK. . . . . . . . . . . . . . . . . . . . . 48 9 CAP´ ITULO 1. INTRODUC¸ ˜ AO 16 1.4 Retornos de ativos Seja Pto pre¸co de um ativo no tempo t. Assumindo que o ativo n˜ao paga dividendos podemos definir os retornos da seguinte forma. Mantendo o ativo num per´ıodo de tempo t´1 at´e o tempo to retorno simples total de um per´ıodo ´e dado por: Rt“Pt Pt´1´1“Pt´Pt´1 Pt´1 (1.13) ou seja 1`Rt“Pt Pt´1 ou Pt“Pt´1p1`Rtq(1.14) Mantendo o ativo em kper´ıodos, entre as datas t´ket, obtemos um retorno simples total de k-per´ıodos 1`Rtrks “ Pt Pt´k˚Pt´1 Pt´2˚¨¨¨˚ Pt´k`1 Pt´k “ p1`Rtqp1`Rt´1q. . . p1`Rt´k`1q “ k´1 ź j“0p1`Rt´jq Na pr´atica o intervalo de tempo ´e muito importante na discuss˜ao dos retornos. Se o intervalo n˜ao for dado explicitamente, ent˜ao assume-se que o mesmo corresponde a um ano. Nesta tese trabalhei com dados di´arios pois se o intervalo for muito grande n˜ao ´e poss´ıvel observar-se alguns fatos estilizados sobre os dados financeiros. Se os ativos forem tomados em kanos, ent˜ao os retornos anuais m´edios s˜ao dados por AnualpRtrksq “ p k´1 ź j“0p1`Rt´jqq1{k´1 (1.15) Outro aspeto importante sobre os retornos ´e a capitaliza¸c˜ao composta dos retornos. A capitaliza¸c˜ao composta dos retornos caracteriza-se por uma fun¸c˜ao exponencial, em que o capital cresce de forma geom´etrica. Muitas vezes por raz˜oes de facilidade de c´alculos ´e usado o logaritmo dos retornos tamb´em conhecido por retorno composto cont´ınuo e ´e dado por: rt“lnp1`Rtq “ ln Pt Pt´1“pt´pt´1(1.16) onde pt“lnpPtq. Uma das vantagens do uso dos logaritmos dos retornos em rela¸c˜ao aos retornos simples ´e que ela permite com que o c´alculo dos logaritmos do retorno dos multi-per´ıodos seja reduzida apenas na soma dos logaritmos dos retornos dos per´ıodos envolvidos. CAP´ ITULO 1. INTRODUC¸ ˜ AO 17 1.5 Carteira de retornos e pagamento de dividendos O retorno simples de uma carteira de Nativos ´e a m´edia ponderada dos retornos simples dos ativos envolvidos, onde o peso em cada ativo ´e a percentagem do valor da carteira investido no ativo. Seja puma carteira que atribui peso wiao ativo i, ent˜ao o retorno simples de pno tempo t ´e Rp,t “ N ÿ i“1 wiRit (1.17) onde Rit ´e retorno simples do ativo i. Seja Dto pagamento de dividendo de um ativo entre as data t´1 e tePto pre¸co do ativo no fim do per´ıodo t. Neste caso Dtn˜ao est´a inclu´ıdo em Pt, os retornos simples e logar´ıtmico s˜ao dados pelas seguintes express˜oes: Rt“Pt`Dt Pt´1´1, rt“lnpPt`Dtq´lnpPt´1q(1.18) Cap´ıtulo 2 Modelos heteroced´asticos Em econometria um dos modelos mais usados para se captar os fatos estilizados apresentados na se¸c˜ao 1.3 s˜ao os modelos ARCH (AutoRegressive Conditional heteroscedasticity) introduzidos por Engle (1982) e generalizados por Bollerslev (1986) GARCH (Generalized Autoregressive Conditional Heteroscedasticity). Estes modelos representam a variˆancia, que ´e o centro da aten¸c˜ao dos modelos, do termo do erro ou inova¸c˜ao atual como quadrado das inova¸c˜oes dos termos anteriores. Esta caracter´ıstica permite que o modelo seja capaz de captar os fatos estilizados dos dados com frequˆencia muito elevadas, por exemplo s´eries temporais financeiras. Neste cap´ıtulo apresentamos esses modelos e as suas varia¸c˜oes mais usadas na pr´atica. 2.1 Modelo GARCH univariado Seja tum processo univariado n˜ao correlacionado com m´edia zero, tsegue um processo ARCH de ordem q, ou simplesmente ARCH(q), se a distribui¸c˜ao condicional de tdado o seu passado Ωt´1“ tt´1, t´2, . . . u, possui m´edia zero e variˆancia condicional dada por σ2 t|t´1“V arpt|Ωt´1q “ V arp2 t|Ωt´1q “ ω`α12 t´1`¨¨¨`αq2 t´p“ω` p ÿ i“1 αi2 t´i(2.1) Na pr´atica, apesar do modelo ARCH captar alguns fatos estilizados das s´eries financeiras como a aglomera¸c˜ao da volatilidade, ´e necess´ario um modelo de ordem muito grande tornando-o pouco parcimonioso. Para colmatar este problema Bollerslev (1986) e Taylor (1986) propuseram um modelo mais geral denominado GARCH com a equa¸c˜ao (2.1) substitu´ıda por σ2 t“ω` q ÿ i“1 αi2 t´i` p ÿ j“1 βjσ2 t´j(2.2) O parˆametro αmede a rea¸c˜ao da volatilidade condicional aos choques do mercado, um valor alto para βindica uma crise nos mercados. A estima¸c˜ao eficiente desses parˆametros ´e uma preocupa¸c˜ao para os profissionais pelo significado que eles desempenham nos mercados de a¸c˜oes. 18 CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 19 Uma abordagem mais abrangente do papel e significado desses parˆametros pode ser encontrada em Carol (2008). Francq et al (2010) afirmam que um processo t´e GARCH(p,q) fraco se este for um ru´ıdo branco (estacion´ario de ordem 4) e o seu quadrado admitir uma representa¸c˜ao ARMA, ou seja 2 t´ r ÿ i“1 ai2 t´i“c`vt´ p ÿ i“1 bivt´i(2.3) onde r = max(p,q) e vt´e a inova¸c˜ao linear de t2 tu. Este resultado ´e muito importante em termos de predi¸c˜ao e sele¸c˜ao de modelos. 2.1.1 Estima¸c˜ao dos parˆametros Existem v´arios m´etodos para a estima¸c˜ao dos parˆametros do modelo GARCH, mas nesta tese apenas se abordou o m´etodo de quase-m´axima verosimilhan¸ca (QML) porque segundo Francq et al (2010) este m´etodo fornece estimadores consistentes e assimptoticamente normais sem exigir a existˆencia de momentos. Seja ttuuma realiza¸c˜ao m,...,nde um processo GARCH(p,q) #εt“σtηt σ2 t“ω`řq i“1αiε2 t´i`řp j“1βjσ2 t´j (2.4) onde ηt´e uma sequˆencia de v.a. iid, na pr´atica ´e comum assumir-se t„Np0,1q,ωą0, αią0 (i“1, . . . , p) e βją0 (j“1, . . . , p). A fun¸c˜ao de QML ´e dada por Lnpθq “ Lnpθ;θ1, . . . , θnq “ n ź t“1 1 a2πσ2 t exp ˆε2 t 2˜σ2 t˙(2.5) onde ˜σt´e obtido recursivamente por (2.2), θ´e um vetor de parˆametros definido por θ“ pθ1, . . . , θp`q`1q“pω, α1, . . . , αq, β1, . . . , βpq. Um estimador QML para θ´e dado por ˆ θn“argmax θPΘ Lnpθq(2.6) onde ˆ θn´e qualquer solu¸c˜ao mensur´avel definida por Θ P p0,`8qxr0,`8qp`q. Maximizar a fun¸c˜ao de m´axima-verosimilhan¸ca ´e equivalente a minimizar, com respeito a θ, ˜ Inpθq “ n´1 n ÿ t“1 ˜ lt onde ˜ lt“˜ ltpθq “ 2 t ˜σ2 t`log ˜σ2 t(2.7) Na pr´atica s˜ao necess´arios valores iniciais para o processo de otimiza¸c˜ao que dependem muito da s´erie em quest˜ao. CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 20 Francq et al (2010) indicam que a variˆancia incondicional dada por varptq “ ω 1´řq i“1αi`řp j“1βj (2.8) ´e um valor razo´avel para se inicializar os valores desconhecidos (2 0,...2 1´qeσ2 0“ ¨¨¨ “ σ2 1´p) Sob certas condi¸c˜oes (para mais detalhes ver Francq et al (2010)) ?npˆ θn´θ0qL ÝÑNp0,pκη´1qJ´1q onde J:“Eθ0ˆd2ltpθ0q dθdθ1˙“Eθ0ˆ1 σ4 tpθ0q dσ2 tpθ0q dθ dσ2 tpθ0q dθ1˙(2.9) ´e uma matriz definida positiva e κn“Erη4 tsă8. Derivando a QML obt´em-se o gradiente dltpθq dθ “ ´ 1 2σ2 t dσ2 t dθ `2 t 2σ4 t dσ2 t dθ “1 2σ2 t dσ2 t dθ ˆ2 t σ2 t´1˙(2.10) com Bσ2 t Bω“1` q ÿ j“1 βjBσ2 t´j Bω(2.11) Bσ2 t Bαi“2 t´i` q ÿ j“1 βjBσ2 t´j Bαpara i“1, . . . , p (2.12) Bσ2 t Bβj“σ2 t´j` q ÿ j“1 βjBσ2 t´j Bβpara j“1, . . . , q et“1, . . . , T. (2.13) Resolver a condi¸c˜ao dσ2 t dθ “0 resulta num conjunto de equa¸c˜oes que podem ser resolvidas pelos m´etodos quase-Newton, por exemplo o algoritmo de Berndt–Hall–Hall–Hausman (BHHH). A itera¸c˜ao BHHH ´e dada por: θt`1“θt`λtH´1 tgt(2.14) onde λt´e um passo vari´avel escolhido para maximizar a QML na dire¸c˜ao apropriada e Ht´e a matriz Hessiana Htpθq “ T ÿ t“1 d2lt dθdθ1(2.15) A matriz hessiana pode ser estimada por Htpθq “ 1 2T T ÿ t“1 ˆ ξtˆ ξ1 t ˆσt 2(2.16) com ξt“ p1, 2 t´1, . . . , 2 t´p, σ2 t´1, . . . , σ2 t´qq. CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 21 Altay et al (2003) prop˜oem uma alternativa dada pelo seguinte problema de optimiza¸c˜ao n˜ao-linear restringido: max ´1 2 T ÿ t“1plog σ2 t`2 t σ2 tq(2.17) ω` q ÿ i“1 αiεt´i` p ÿ j“1 βjσ2 t´j“σ2 t@t“1, . . . , T m ÿ i φiYt´i`t“Yt@t“1, . . . , T σ2 tě0@t“1, . . . , T αiě0@i“1, . . . , q βjě0@j“1, . . . , p ωą0 a qual denominam de modelo GARCH-NLP, onde řm iφiYt´i`t“Yt´e a equa¸c˜ao autorregressiva usada para representar os dados. Segundo Andersen (2009) as estimativas encontradas sob o uso da QML s˜ao assint´oticamente normalmente distribu´ıdas desde que a m´edia condicional e a fun¸c˜ao da variˆancia condicional do modelo GARCH sejam bem especificadas. Bollerslev e Wooldridge (1992) resolveram o problema do c´alculo dos erros padr˜ao das estimativas derivando uma matriz de covariˆancia robusta `a n˜ao normalidade para QML dada por Hpˆ θQMLq´1Bpˆ θQMLqHpˆ θQMLq´1(2.18) onde ˆ θQML denota a QMLE de θe B ´e dada por Bpθq “ T ÿ t“1 Blt BθBlt Bθ1.(2.19) Esta foi a matriz de covariˆancia usada para se obterem os erros padr˜ao apresentados nesta tese. 2.1.2 Testando os efeitos ARCH Antes de se prosseguir para a estima¸c˜ao ou sele¸c˜ao de um modelo ARCH ´e pr´atica comum ser conduzido um teste sobre a existˆencia de efeitos ARCH ou heteroced´asticos nos res´ıduos. O teste ´e usualmente feito nos quadrados da s´erie. Para conduzir esse teste temos duas possibilidades. A primeira consiste em aplicar a estat´ıstica Qpmqde Ljung-Box nos quadrados da s´erie, sendo a hip´otese nula de que a s´erie ´e n˜ao correlacionada nos seus primeiros mdesfasamentos. A outra possibilidade ´e usar o teste multiplicativo de Lagrange proposto por Engle. Basicamente este teste ´e similar ao teste Fconduzido `a regress˜ao linear CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 22 X2 t“α0`α1X2 t´1`¨¨¨`αmX2 t´m`t, T “m`1, . . . , T sendo a estat´ıstica F dada por F“pSSR0´SSR1q{m SSR1{pT´2m´1q(2.20) que sob a hip´otese nula possui uma distribui¸c˜ao assintoticamente qui-quadro com mgraus de liberdade, onde SSR0“ T ÿ t“m`1pX2 t´¯ωq2,¯ω“1 T T ÿ t“1 2 teSSR1“ÿ t“1`m Tˆ2 t. Rejeita-se a hip´otese nula se Fąχ2 αse o păαsendo χ2 mpαqo percentil 100p1´αqdo quantil superior de χ2 m. Outros modelos foram propostos na literatura para se captar a dinˆamica das s´eries com elevada frequˆencia como ´e caso de vari´aveis macroecon´omicas e de alguns fatos estilizados das s´eries financeiras. A seguir apresentamos algumas das especifica¸c˜oes mais usadas na pr´atica. 2.1.3 Outros modelos GARCH Modelo GJR-GARCH Este modelo foi proposto por Glosten et al (1993). Neste modelo introduz-se um novo parˆametro para captar o efeito de ”leverage”quando o choque nos mercados s˜ao negativos: σ2 t“ω` q ÿ i“1 αiεt´i`λIrt´iă0s` p ÿ j“1 βjσ2 t´j(2.21) onde a fun¸c˜ao indicadora Irt´iă0s“1 se tă0 e 0 caso contr´ario. A estima¸c˜ao dos parˆametros ´e feita usando a fun¸c˜ao de QML (2.5). Modelo EGARCH Este modelo garante que a variˆancia seja positiva sem impor restri¸c˜oes nos coeficientes, formulando a equa¸c˜ao da variˆancia condicional em termos do seu logaritmo. Para um processo EGARCH(p,q) a variˆancia condicional ´e modelada pela seguinte express˜ao: log σ2 t“ω` q ÿ i“1 αigpZt´iq` p ÿ j“1 βjlog σ2 t´j(2.22) onde gpZtq “ θZt`λp|Zt|´Ep|Zt|q,ω, α, β s˜ao constantes definidas como anteriormente e Ztpode ser uma vari´avel normal padr˜ao ou uma distribui¸c˜ao de erro generalizada. CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 23 A fun¸c˜ao de m´axima-verosimilhan¸ca para este modelo ´e dada por: ln Lpω, θ, λ, βq“´1 2 T ÿ t“1plnpσ2 tq`pt σtq2(2.23) 2.2 Modelos GARCH multivariados Carol (2008) defende que para al´em da aglomera¸c˜ao da volatilidade resultante das m´as not´ıcias que podem surgir nos mercados, tamb´em existe aquela que ´e causada pela correla¸c˜ao dos pre¸cos de ativos. Para capt´a-la ´e necess´ario estimarmos a matriz de covariˆancia, usando o modelo GARCH multivariado. Existem v´arias extens˜oes do modelo GARCH multivariado. Abaixo apresentamos os dois modelos mais gerais e dos primeiros a serem formulados, uma vez que muitas outras especifica¸c˜oes s˜ao derivadas deles. 2.2.1 Modelo VECH O modelo VECH ´e a formula¸c˜ao multivariada direta do modelo univariado. Para Francq et al (2010) um processo ttu“t1,...,Nu´e um VEC-GARCH(p, q) quando satisfazer #εt“H 1 2 tηt vechpHtq “ ω`řq i“1Apiqvechpεt´iε1 t´iq`řp j“1BpjqvechpHt´jq(2.24) onde vech (.) ´e o operador que consiste em empilhar os elementos das por¸c˜oes da matriz triangular inferior e diagonais de uma matriz sim´etrica num vetor coluna, Ht´e uma matriz definida positiva, ηt´e uma sequˆencia de v.a iid com distribui¸c˜ao η,ω´e um vetor de tamanho tNpN`1q{2ux 1 , ApiqeBpjqs˜ao matrizes de dimens˜ao {NpN`1q{2}x{NpN`1q{2}. Este modelo possui grandes limita¸c˜oes na pr´atica, porque n˜ao ´e poss´ıvel derivar diretamente a propriedade definida positiva de Hte ´e muito custoso a n´ıvel computacional uma vez que mesmo para ordens pequenas ´e necess´ario estimar um elevado n´umero de parˆametros. Se ttusegue uma distribui¸c˜ao normal a fun¸c˜ao de log-m´axima-verosimilhan¸ca do modelo GARCH multivariado ´e dada por T ÿ t“1 lt“1 2Nlnp2πq´ 1 2ln |Ht|´ 1 2 T ÿ t“1 1 tH´1 t´1t(2.25) Para facilitar o tratamento deste modelo s˜ao impostas algumas condi¸c˜oes nas matrizes Apiq eBpjq. O modelo proposto nesta tese assume que as matrizes ApiqeBpjqs˜ao sim´etricas. A seguir apresentamos as especifica¸c˜oes mais comuns do modelo VECH que s˜ao os modelos DVECH e BEEK. CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 24 2.2.2 Modelo DVECH O modelo DVECH (Diagonal VECH) permite que uma representa¸c˜ao vetorial seja obtida pelas matrizes diagonais ApiqeBpjq. A variˆancia condicional ´e dada por hkl,t “ωkl ` q ÿ i“1 apiq kl k,l´ik,l´i` p ÿ j“1 bpjq kl hkl,t´j(2.26) com ωkl “ωlk,apiq kl “apiq lk ebpjq kl “bpjq lk @pl, kq. Para N“1, este modelo ´e equivalente ao modelo GARCH univariado, para o caso em que Ną1, a equa¸c˜ao (2.26) pode ser escrita da seguinte forma matricial: Ht“¨ ˚ ˝ ω11 ¨¨¨ ω1m . . .. . . ω1m¨¨¨ ωmm ˛ ‹ ‚` q ÿ i“1¨ ˚ ˝ apiq 11 2 1,t´i¨¨¨ apiq 1N1,t´iN,t´i . . .. . . apiq 1m1,t´im,t´i¨¨¨ apiq NN 2 m,t´i ˛ ‹ ‚(2.27) ` p ÿ j“1¨ ˚ ˝ bpjq 11 h1,t´i¨¨¨ bpjq 1Nh1N,t´i . . .. . . bpjq 1mh1,t´i¨¨¨ bpjq NN 2 N,t´i ˛ ‹ ‚ :“Ω` q ÿ i“1 diagpt´iqApiqdiagpt´iqApiq` p ÿ j BpjqdHt´j onde d´e o produto de Hadamard, isto ´e, elemento por elemento. 2.2.3 Modelo BEEK O modelo BEEK ´e uma especifica¸c˜ao do modelo multivariado proposto por Baba, Engle, Kraft e Kroner. Segundo Francq et al (2010) um modelo GARCH-BEEK(p,q) ´e todo o processo ttucom respeito a uma sequˆencia de v.a. iid ηtque satisfaz #εt“H 1 2 tηt Ht“Ω`řq i“1řK k“1Aikεt´iε1 t´iA1 ik `řp j“1řK k“1BjkHt´jB1 jk (2.28) onde K´e um inteiro, Ω, Aik eBjk s˜ao matrizes NxNe Ω ´e definida positiva. 2.2.4 Estacionariedade do modelo VECH A estacionariedade do modelo VECH ´e pouco consensual dentro da comunidade cient´ıfica. Na maior parte da literatura sobre os modelos VECH a estacionariedade ´e obtida por meio da representa¸c˜ao ARMA do modelo. Seja a dinˆamica da volatilidade do modelo dada por (2.24) considerando vt“vechpttq ´ vechpHtqum processo de ru´ıdo branco e reorganizando os termos em (2.24) obt´em-se um processo autogressivo de m´edia m´ovel (VARMA) de ordem (maxpp, qq, q) escrevendo o processo de ordem 2 ηtda seguinte forma CAP´ ITULO 2. MODELOS HETEROCED ´ ASTICOS 25 ηt“ω` maxpp,qq ÿ i“1pApiq`Bpiqqvechpttq`vt´ p ÿ j“1 Bpjqvt´j(2.29) onde Ap`1“ ¨¨¨ “ Aq“0 se pąqeBq`1“ ¨¨¨ “ Bp“0 se qąp. A partir das caracter´ısticas do modelo ARMA pode-se verificar que o modelo (2.24) ´e estacion´ario de segunda ordem se, e somente se, os valores pr´oprios em m´odulo de maxpp,qq ÿ i“1 Apiq`Bpiq(2.30) s˜ao menores que 1. Francq et al (2010) s˜ao de opini˜ao que o resultado apresentado acima n˜ao ´e de todo correto uma vez que vechpttq “ ˜r ÿ i“1pAi`BiqBi¸´1˜ω`vt´ p ÿ j“1 Bjvt´j¸(2.31) apenas constitui uma solu¸c˜ao se vt“vechpttq ´ vechpHtqpoder ser expressa como uma fun¸c˜ao de {ηt´u, u ą0}. Apesar deste resultado referido por Francq et al (2010), na pr´atica verificou-se que quando as matrizes ApiqeBpiqobedeciam (2.30) os modelos obtidos a partir da simula¸c˜ao dessas matrizes eram estacion´arios. O tratamento das condi¸c˜oes de estacionariedade para o modelo BEKK e DVECH podem ser consultadas em Francq et al (2010). 2.2.5 Estima¸c˜ao dos parˆametros A express˜ao abaixo ´e a extens˜ao multivariada da alternativa GARCH-NLP proposta por Altay (2003): max ´1 2 T ÿ t“1plog detHt`T tH´1 ttq vechpHtq “ ω` q ÿ i“1 Apiqvechpεt´iε1 t´iq` p ÿ j“1 BpjqvechpHt´jq @t“1, . . . , T m ÿ i φiYl,t´i`lt “Ylt l“1, . . . , l @t“1, . . . , T Htľ0@t“1, . . . , T onde o s´ımbolo ľsignifica ”sim´etrica, definida positiva”. CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 32 0 200 400 600 800 1000 −60 −20 20 60 Volatilidade anualizada de SP500 n volatiliade 0 200 400 600 800 1000 −50 0 50 Volatilidade anualizada de FTSE n volatiliade Figura 3.6: Volatilidade condicional anualizada vs retornos anualizados das ´ultimas 500 observa¸c˜oes. Na figura 3.6 podemos verificar que os retornos do FTSE s˜ao mais vol´ateis em rela¸c˜ao aos retornos do S&P500. Apesar do modelo GARCH(1,1) sim´etrico ser o mais f´acil de se aplicar, e na realidade o mais usado na pr´atica, ele n˜ao consegue captar o efeito de leverage como visto na figura 3.6, uma vez que ela n˜ao leva em conta os valores negativos dos retornos. Esses resultados comparados com outros pacotes encontrados nos software R e Matlab mostraram muita similitude uma vez que todos eles usam a mesma metodologia. A raz˜ao da ligeira diferen¸ca ´e devido a certos parˆametros iniciais, a implementa¸c˜ao e a diferentes rotinas de optimiza¸c˜ao usadas na estima¸c˜ao dos parˆametros. O caso univariado foi usado apenas como ensaio e para obter uma melhor percep¸c˜ao da dinˆamica de cada carteira com finalidade de estudar a intera¸c˜ao da mesma de forma conjunta. Os coeficientes βS&P500 “0.961793 e βF T SE “0.902044 indicam que a variˆancia condicional leva muito tempo a extinguir-se, ou seja, a volatilidade ´e persistente, ao passo que os coeficientes αS&P500 “0.031228 e αF T SE “0.063935 indicam que a volatilidade possui uma rea¸c˜ao pouco intensiva aos movimentos dos mercados. 3.3 Modela¸c˜ao da volatilidade dos retornos de ´ındices atrav´es dos modelos GARCH: caso multivariado Para o caso multivariado usou-se o mesmo conjunto de dados onde a janela temporal foi escolhida usando os teste de Engle sobre a existˆencia de efeitos ARCH. Foi aplicado o modelo VECH apenas com a restri¸c˜ao de que as matrizes ApiqeBpjqsejam sim´etricas usando a formula¸c˜ao GARCH NLP, ou mais simplesmente NLP, proposta por Altay et al (2003). CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 33 Os c´alculos foram realizados com o pacote fmincon da toolbox de optimiza¸c˜ao do Matlab vers˜ao 8.1 sendo os modelos DVEC e BEEK usados como modelos competidores sendo estimados com o pacote FINMETRICS 2.0 do S-Plus vers˜ao 8. De salientar que quanto maior for o n´umero de parˆametros, mais achatada ´e a superf´ıcie da fun¸c˜ao de m´axima verosimilhan¸ca (Carol, 2008) o que conduz `a existˆencia de v´arios pontos m´ınimos/m´aximos locais. Para minimizar este problema foi tamb´em usado o pacote de optimiza¸c˜ao global do Matlab para efetuar uma busca a partir de v´arios pontos. Um dos problemas recorrentes na estima¸c˜ao do modelo multivariado ´e encontrar ”bons” parˆametros iniciais que assegurem que o procedimento tende assintoticamente para os parˆametros ´optimos. Uma solu¸c˜ao ”ad hoc”´e A“ p?0.15 ´cqI`cˆJ B“ p?0.8´cqI`cˆJ, onde I´e uma matriz identidade da ordem de A,J´e uma matriz quadrada da ordem de A com elementos todos iguais a 1 e cum n´umero real muito pequeno, por exemplo 0.05. Os valores 0.15 e 0.8 foram escolhidos uma vez que no caso GARCH(1,1) univariado os mesmos serem os principais candidatos para se inicializar os parˆametros αeωrespetivamente. Nota-se que esses valores devem garantir que o processo seja estacion´ario, isto ´e, a sua soma deve ser menor que 1, de salientar que outros valores poderiam ser usados desde que os mesmos gerem um processo estacion´ario GARCH(1,1) univariado. A matriz C´e inicializada com a decomposi¸c˜ao Cholesky de 20(“1 1´0.15´0.8) vezes a covariˆancia n˜ao condicional dos dados. O problema da inicializa¸c˜ao da matriz de variˆancia condicional Htfoi resolvido inicializandoa com a variˆancia n˜ao condicional dos dados. CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 34 Parˆametros DVEC BEKK NLP C(1,1) 0.006101 0.057854* 0.0098** (0.016081) (0.008720) (0.0046) C(2,1) 0.013974 0.046324 0.0126 (0.020129) (0.060460) (0.0076) C(2,2) 0.02306* 0.174900* 0.0234* (0.007302) (0.026505) (0.0130) A(1,1) 0.034965* 0.153829* 0.0288** (0.006137) (0.015968) (0.0091) A(2,1) 0.019760 -0.053035* 0.0109 (0.011796) (0.024387) (0.0063) A(1,2) 0.018800 (0.020850) A(1,3) 0.0116** (0.0071) A(2,2) 0.051744* 0.254250* 0.0066 (0.009256) (0.028261) (0.0143) A(2,3) 0.0226* (0.0102) A(3,3) 0.0590** (0.0223) B(1,1) 0.952451* 0.984434* 0.9589** (0.00753) (0.003495) (0.0124) B(2,1) 0.890176* 0.012521 -0.0049** (0.081961) (0.007105) (0.0102) B(1,2) -0.002573 (0.008268) B(1,3) -0.0189* (0.0092) B(2,2) 0.911816* 0.940532* 0.8958 (0.018914) (0.013902) (0.0564) B(2,3) -0.0187** (0.0193) B(3,3) 0.9074 (0.0324) Verossimilhan¸ca -3205.814 -3215.106 -3195.8 AIC 6433.629 6456.211 6421.6 BIC 6492.015 6525.214 6501.3 Tabela 3.4: Resultado do caso bivariado sobre retornos de ´ındices S&P500 e FTSE (valores entre parˆenteses representam os erros padr˜ao), (*) coeficientes significativos com erros padr˜ao robustos e (**) coeficientes significativos com erros padr˜ao robustos e obtidos por aproxima¸c˜ao num´erica. CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 35 Tendo em conta os parˆametros significativos na ´ultima coluna da tabela 3.4 , os elementos da matriz de covariˆancia s˜ao representados pelas seguintes equa¸c˜oes recursivas: hS&P500,t “0.0098 `0.02882 S&P500,t´1`0.01162 F T SE,t´1`0.9589hS&P500,t´1´0.01189hF T SE,t´1 hS&P500ˆF T SE,t “0.02262 F T SE,t´1`0.8958hS&P500ˆF T SE,t´1 hF T SE,t “0.0226S&P500,t´1F T SE,t´1`0.05902 F T SE,t´1´0.0189hS&P500,t´1`0.9109hF T SE,t´1 (3.2) Verifica-se claramente que a volatilidade cruzada entre os retornos ´e apenas modelada pelos quadrados dos retornos FTSE e por si mesma, o que como veremos a seguir n˜ao permite que o modelo capte a volatilidade cruzada dos retornos em quest˜ao. A tabela 3.4 mostra os resultados da estima¸c˜ao dos modelos DVECH, BEKK e NLP com os seus respectivos coeficientes estimados e os valores dos crit´erio de informa¸c˜ao. Os crit´erios de informa¸c˜ao e a fun¸c˜ao de quase m´axima-verosimilhan¸ca mostram que o modelo NLP possui melhor desempenho em rela¸c˜ao aos modelos competidores. Figura 3.7: Volatilidade de S&P500 Nas figuras 3.7e3.8 apresentam-se os gr´aficos da volatilidade anualizada de S&P500 e FTSE obtidos pelos trˆes m´etodos multivariados. As express˜oes usadas s˜ao semelhantes `as usadas no caso univariado. CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 36 Figura 3.8: Volatilidade de FTSE Quanto mais pr´oximo o modelo estiver das linhas descont´ınuas melhor ´e o ajustamento aos dados e melhor consegue explicar e captar a dinˆamica do processo. A visualiza¸c˜ao desses gr´aficos permite-nos verificar que os modelos DVECH e BEKK permanecem quase constantes n˜ao captando assim a volatilidade dos retornos, ao passo que o modelo NLP (restringida) parece adequar-se melhor ao padr˜ao da dinˆamica dos retornos no per´ıodo de tempo observado. O gr´afico da covariˆancia condicional entre os ´ındices S&P500 e FTSE foi obtido multiplicando os retornos de S&P500 pelos retornos de FTSE. O resultado ´e apresentado na figura 3.9, nela podemos verificar a dificuldade dos modelos em captar a covariˆancia ou dependˆencia correlacional entre os ´ındices. Contudo o modelo NLP tamb´em possui um desempenho ligeiramente melhor em rela¸c˜ao aos outros modelos competidores por causa das demasiadas condi¸c˜oes que neles s˜ao impostas. 3.3.1 Avalia¸c˜ao do modelo Carol (2008) afirma que um modelo GARCH captura a volatilidade efetivamente quando os retornos padronizados pela sua variˆancia condicional possuem uma distribui¸c˜ao normal. E, para a mesma autora, o modelo ´e bem especificado se o quadrado dos res´ıduos padronizado pela variˆancia condicional n˜ao possuir correla¸c˜oes. Essas ideias foram postas em pr´atica abaixo para se testar a qualidade do ajustamento do modelo NLP. A avalia¸c˜ao dos modelos competidores BEKK e DVECH podem ser vistas no apˆendice B. Os correlogramas dos res´ıduos padr˜ao (figura 3.10) obtidos pelo modelo NLP mostra que CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 37 Figura 3.9: Covariˆancia condicional de S&P500 e FTSE Figura 3.10: Correlograma dos res´ıduos padronizados este modelo se ajusta `a estrutura correlacional da s´erie no modelo condicional, uma vez que o correlograma indica n˜ao haver correla¸c˜ao nas s´eries ap´os o ajustamento do modelo aos CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 38 dados. SP&500 lags h p-valor estat´ıstica valor cr´ıtico 10 0 0.9040 4.8015 18.3070 15 0 0.7849 10.5343 24.9958 20 0 0.9021 12.3876 31.4104 FTSE lags h p-valor estat´ıstica valor cr´ıtico 10 0 0.9377 4.2032 18.3070 15 0 0.8626 9.2761 24.9958 20 0 0.8404 13.8021 31.4104 Tabela 3.5: Resultado do teste sobre a n˜ao correla¸c˜ao dos res´ıduos padronizados de S&P500 e FTSE. Para formalizar as conclus˜oes conduziu-se o teste de Ljung-Box cuja hip´otese nula ´e que a s´erie ´e n˜ao correlacionada, sobre o quadrado dos res´ıduos padronizados. A tabela 3.5 mostra os resultados deste teste aplicado ao quadrado dos res´ıduos padronizados, nela podemos verificar o n´umero de desfasagem, o p-valor, a estat´ıstica do teste e o valor cr´ıtico. O hrepresenta uma vari´avel indicadora bin´aria onde o valor 0 indica n˜ao rejei¸c˜ao da hip´otese nula e 1 o caso contr´ario. Com base nos resultados apresentados nessas tabelas podemos concluir formalmente (quantitativamente) que os quadrados dos res´ıduos padronizados s˜ao n˜ao correlacionados, por conseguinte o modelo NLP captura a aglomera¸c˜ao da volatilidade nos retornos dos ´ındices S&P500 e FTSE. p-valor FValor cr´ıtico S&P500 0.8182 0.0529 3.8415 FTSE 0.5480 1.2030 5.9915 Tabela 3.6: Teste dos efeitos ARCH sobre os retornos de padronizados. Os p-valor na tabela 3.6 mostra que n˜ao h´a evidˆencias estat´ısticas para se rejeitar a hip´otese nula de que n˜ao haja efeitos ARCH sobre os retornos padronizados sob a hip´otese alternativa de um modelo ARCH(2). CAP´ ITULO 3. APLICAC¸ ˜ AO DOS MODELOS GARCH 39 Figura 3.11: QQ plot dos res´ıduos padronizados O QQ-plot representado na figura 3.11 mostra que os res´ıduos padronizados se afastam da normalidade, principalmente nas caudas, algo que ´e confirmado pelos elevados valores da estat´ıstica de Jarque–Bera na tabela 3.7 rejeitando-se a hip´otese nula de que os res´ıduos provˆem de uma distribui¸c˜ao normal. p-valor FValor cr´ıtico S&P500 0.0010 253.4749 5.9513 FTSE 0.0010 137.6875 5.9513 Tabela 3.7: Teste sobre a normalidade dos retornos padronizados. Estes resultados mostram claramente que, tal como no caso univariado, muitas vezes ´e necess´ario uma distribui¸c˜ao com caudas mais pesadas como ´e o caso da distribui¸c˜ao t sim´etrica. Carol (2008) afirma que quando dois modelos GARCH se ajustam satisfatoriamente aos dados, escolhe-se o que gerar maior valor para fun¸c˜ao de m´axima-verosimilhan¸ca. Tendo em conta isso, podemos concluir que o NLP ´e o modelo que melhor modela a volatilidade dos retornos dos ´ındice S&P500 e FTSE. Cap´ıtulo 4 Conclus˜oes A computa¸c˜ao do modelo NLP no Matlab usando as ferramentas de acelera¸c˜ao dos scripts (que envolveu a convers˜ao do c´odigo para linguagem C++) teve um tempo de execu¸c˜ao em m´edia n˜ao superior a 3 minutos, ao passo que usando os algoritmos de optimiza¸c˜ao global , com o mesmo resultado, o tempo de execu¸c˜ao foi praticamente 5 vezes superior em rela¸c˜ao aos algoritmos de optimiza¸c˜ao local uma vez que a busca do ponto ´otimo ´e feita a partir de v´arios pontos iniciais selecionando-se o candidato que gerar o menor dos pontos m´ınimo. No S-Plus a estima¸c˜ao dos modelos DVECH e BEKK levou apenas 1/3 do tempo observado nos scripts do Matlab. O comportamento do ´ındice FTSE afeta significativamente o comportamento do ´ındice S&P500 quer em termos de retornos quer em termos de volatilidade. Ganhos de FTSE afetam favoravelmente o S&P500. Aumento ou diminui¸c˜ao da volatilidade de FTSE afeta no sentido contr´ario o S&P500. A “persistˆencia” da volatilidade ´e mais alta no ´ındice S&P500 do que no FTSE. Quanto maior for a persistˆencia na volatilidade mais lentamente o efeito de choques externos na volatilidade s˜ao dissipados no tempo. Apesar disso, nota-se que a volatilidade de FTSE atinge em v´arios momentos valores mais altos do que os do S&P500. O modelo falha em captar a correla¸c˜ao cruzada dos retornos, a inclus˜ao do quadrado dos retornos de FTSE n˜ao melhorou o modelo. O mesmo aconteceu com os modelos competidores BEKK e DVECH por causa da pouca significˆancia dos parˆametros que modelam a volatilidade cruzada dos retornos dos ´ındices S&P500 e FTSE. Os resultados mostram que o modelo NLP restringido ´e uma alternativa aos modelos tradicionais DVECH e BEKK a ter em conta na solu¸c˜ao dos problemas de estima¸c˜ao do modelo GARCH, tendo representado melhor a volatilidade dos ´ındice S&P500 e FTSE para o intervalo de tempo usado que os modelos BEKK e DVECH. 40 Bibliografia 41 APˆ ENDICE B. RESULTADOS DOS MODELOS DVECH E BEKK. 48 Figura B.2: Correlograma dos res´ıduos padronizados do modelo DVECH F-estat´ıstica p-valor chi2df S&P500 6.612 0.8822 12 FTSE 5.826 0.9246 12 Tabela B.5: Teste de Ljung-Box sobre a existˆencia de autocorrela¸c˜ao nos retornos padronizados do modelo BEKK. estat´ıstica p-valor S&P500 6.6086 0.8824 FTSE 5.8198 0.9249 Tabela B.6: Teste de multiplicador de Lagrange sobre existˆencia de efeitos ARCH nos retornos padronizados do modelo BEKK. APˆ ENDICE B. RESULTADOS DOS MODELOS DVECH E BEKK. 49 Figura B.3: QQ plot dos res´ıduos padronizados do modelo BEKK Figura B.4: Correlograma dos res´ıduos padronizados do modelo BEKK Apˆendice C C´odigos Matlab 50 APˆ ENDICE C. C ´ ODIGOS MATLAB 51 1function [parametros, ll, Ht, VCV, scores,HESSIAN] = mgvech(parametros iniciais,... 2dados, p, q,tipo,F) 3%MGVECH estima o modelo VECH-GARCH 4% 5% Inputs: 6% parametros iniciais: uma matriz coluna representando os parametros 7% iniciais do modelo dados por $$c 1, c 2, c 3, a {1,1}, a {2,2}, 8% a {3,3}..., a {n,n}, , a {1,2},, a {1,3}..., a {m,n},b {1,1}, 9% b {2,2}, b {3,3}..., b {n,n}, , b {1,2},, b {1,3}..., b {m,n}$$ 10 % 11 % dados: matriz k x T T $$ de v.a iid com zero e variancia *I*. o 12 %*k*representa o numero de retornos de ativos e T numero de 13 % observacoes 14 % 15 % p: ordem das inovacoes ou o numero de defasagem dos coeficientes ARCH. 16 % 17 % q: ordem das variancias condicionais ou o numero de defasagem dos 18 % coeficientes GARCH. 19 % 20 % sinal: usada tornar negativa a funcao de maximo-verosimilhanca default 21 % -1 22 % 23 % F: (function handle) usada no caso de acelerarmos as simulacoes usando 24 % um ficheiro mex. 25 % 26 % Outputs: 27 % 28 % 29 % parametros = MGVECH(parametros, dados, p, q, sinal) estima os 30 % parametros (*AeB*) do modelo. 31 % 32 % [parametros, ll] = MGVECH(parametros, dados, p, q, sinal,F) 33 % retorna o valor da maximo verosimilhanca avaliada nos *parametros* 34 % 35 % [parametros, ll, Ht, VCV] = MGVECH(parametros, dados, p, q, sinal,F) 36 % retornas a matriz das variancias condicionais (volatilidade). 37 % 38 % [parametros, ll, Ht, VCV] = MGVECH(parametros, dados, p, q, sinal,F) 39 % retornas a matriz das covariancias condicionais. 40 % 41 % [parametros, ll, Ht, VCV, scores] = MGVECH(parametros, dados, p, q, 42 % sinal,F) retornas a matriz dos scores. 43 44 45 [dadosSim, backcast,T,k] = data transform(dados); 46 47 %limites dos parametros 48 %dims = k*(k+1)/2; 49 UB = .99998 *ones(size(parametros iniciais)); 50 LB = -.99998 *ones(size(parametros iniciais)); 51 52 %LB(1:dims) = eps; 53 %LB(dims+1) = eps; 54 %LB(dims*dims+1) = eps; 55 APˆ ENDICE C. C ´ ODIGOS MATLAB 52 56 cons = @constraint; 57 if nargin <6 58 F = @likelihoods; 59 end 60 61 parametros = []; 62 HESSIAN = []; 63 64 if strcmpi(tipo,'normal') 65 options = optimoptions('fmincon','Algorithm','interior-point',... 66 'MaxFunEvals',6000); 67 options = optimoptions(options,'display','iter'); 68 [parametros,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN] = fmincon(F,parametros iniciais,[],[],[],[],LB,UB,cons,... 69 options,dados,dadosSim,backcast,p,q,T,k); 70 HESSIAN = HESSIANˆ(-1); 71 warning('on') 72 73 elseif strcmpi(tipo,'global') 74 F = @(parameters)F(parameters,dados,dadosSim,backcast,p,q,T,k); 75 cons = @(parameters)constraint(parameters,dados,dadosSim,backcast,... 76 p,q,T,k); 77 LB = eps *ones(size(parametros iniciais)); 78 % Set options for FMINCON 79 options = optimset('Algorithm','interior-point','Disp','iter',... 80 'MaxFunEvals',6000); 81 % Create problem for GlobalSearch 82 problem = createOptimProblem('fmincon','objective',F,'x0',... 83 parametros iniciais,'lb',LB,'ub',UB,'nonlcon',cons,... 84 'options',options); 85 86 % Make a GlobalSearch object 87 gs = GlobalSearch; 88 89 % Run the optimization 90 [parametros,fval,eflag,output,solutions] = run(gs, problem); 91 warning('on') 92 end 93 94 %passar o oposto de sinal uma vez que precisamos o valor que maximiza 95 %a funcao. 96 F = @likelihoods; 97 [ll,~,Ht] = F(parametros, dados, dadosSim, backcast, p, q, T, k); 98 ll = -ll; 99 100 %Realizar a inferencia do modelo 101 102 if nargout>=4 103 [VCV,~,~,scores] = robustvcv(F,parametros,0,dados, dadosSim, ... 104 backcast,p,q,T,k); 105 end 106 107 end APˆ ENDICE C. C ´ ODIGOS MATLAB 53 1function [lik, liks, H] = likelihoods(parameters, data, datasym, backcast,p,q,T,k) 2%LIKELIHOODS calcula a funcao de maximo-verosimilhanca do modelo garch 3% multivariado. 4% 5% Inputs: 6% parametros: uma matriz coluna representando os parametros 7% iniciais do modelo dados por $$c 1, c 2, c 3, a {1,1}, a {2,2}, 8% a {3,3}..., a {n,n}, , a {1,2},, a {1,3}..., a {m,n},b {1,1}, 9% b {2,2}, b {3,3}..., b {n,n}, , b {1,2},, b {1,3}..., b {m,n}$$ 10 % 11 % dados: matriz k x T T $$ de v.a iid com zero e variancia *I*. o 12 %*k*representa o numero de retornos de ativos e T numero de 13 % observacoes. 14 % 15 % datasym: matriz k x k x T de v.a iid com zero e variancia 16 %*I*. o *k*representa o numero de retornos de ativos e T numero 17 % de observacoes. 18 % 19 % backcast: matriz k x T T $$ usada para inicializar a variancia 20 % condicional no caso do numero de observacao ser menor que p ou q, o 21 % valor pre-definido e covariancia nao condicional dos dados. 22 % 23 % p: ordem das inovacoes ou o numero de defasagem dos coeficientes ARCH. 24 % 25 % q: ordem das variancias condicionais ou o numero de defasagem dos 26 % coeficientes GARCH. 27 % 28 % sinal: usada tornar negativa a funcao de maximo-verosimilhanca default 29 % -1 30 % 31 % Outputs: 32 % 33 % lik = likelihoods(parametros, dados, dadosSim, backcast,p,q,T,k,sinal) 34 % retorna o valor da funcao de maximo-verosimilhanca avaliados no para35 % mentros. 36 % 37 % [lik, liks] = likelihoods(parametros, dados, dadosSim,backcast,p,q,T, 38 % k,sinal) retorna o vetor das maximo-verosimilhanca obtida em cada 39 % observacao (ou itera??o). 40 % 41 % [lik, liks, H] = likelihoods(parametros, dados, dadosSim, backcast,p,q, 42 % T,k,sinal) retorna a matriz das variancias condicionadas(volatilidade). 43 % 44 % Referencias 45 % Christian Francq, Jean-Michel Zakoian (2010), GARCH Models: Structure, 46 % Statistical Inference and Financial Applications,Wiley. 47 % 48 % Torben G. Andersen, Richard A. Davis Jens-Peter Kreig,Thomas Mikosch 49 % (2009), Handbook of Financial Time Series. 50 % 51 % Brooks, Chris (2008) Introductory Econometrics for Finance 2nd EDITION. 52 53 54 if size(data,1) >size(data,2) 55 data = data'; APˆ ENDICE C. C ´ ODIGOS MATLAB 54 56 end 57 58 [C,A,B] = par transform(parameters,k); 59 60 H = zeros(k,k,T); 61 C = unvech(C); 62 liks = zeros(T,1); 63 64 65 c = .5*k*log(2*pi); 66 67 for t = 1:T 68 H(:,:,t) = C; 69 for i = 1:p 70 if t<= i 71 H(:,:,t) = H(:,:,t) + unvech(A(:,:,i)*vech(backcast)); 72 else 73 H(:,:,t) = H(:,:,t) + unvech(A(:,:,i)*vech(datasym(:,:,t-i))); 74 end 75 end 76 77 for i = 1:q 78 if t<= i 79 H(:,:,t) = H(:,:,t) + unvech(B(:,:,i)*vech(backcast)); 80 else 81 H(:,:,t) = H(:,:,t) + unvech(B(:,:,i)*vech(H(:,:,t-i))); 82 end 83 end 84 %the line bellow is just used to make a mex function 85 aux = c +.5 *log(complex((det(H(:,:,t))))) + 0.5*data(:,t)'*(H(:,:,t)ˆ(-1))*data(:,t); 86 if imag(aux) == 0 87 liks(t) = real(aux); 88 else 89 liks(t) = nan; 90 end 91 % liks(t) = c - .5 *log(det(H(:,:,t))) - 0.5*data(:,t)'*(H(:,:,t)ˆ(-1))*data(:,t)); 92 end 93 94 lik = sum(liks); 95 if isnan(lik) | | isinf(lik) | | ~isreal(lik) 96 lik = 1e7; 97 end APˆ ENDICE C. C ´ ODIGOS MATLAB 55 1function [c,ceq] = constraint(parameters, data, datasym, backcast,p,q,T,k) 2 3 4[~, A, B] = par transform(parameters,k); 5 6c = abs(eig(A+B)) - 1 + eps; 7ceq = []; APˆ ENDICE C. C ´ ODIGOS MATLAB 56 1function [C, A, B, par] = generate par(data,a,b,c,p,q,k) 2 3 4if nargin <2 5c = 0.05; 6p=1; q=1; 7a = .15; 8b = .8; 9k = 2; 10 end 11 12 d=k*(k+1)/2; 13 14 I = eye(d); 15 J = ones(d); 16 17 A = zeros(d,d,p); 18 B = zeros(d,d,q); 19 20 scale = 100; 21 for i=1:p 22 A(:,:,i) = ((sqrt(a) - c)*I+c*J)/scale; 23 end 24 25 for j=1:q 26 B(:,:,j) = ((sqrt(b) - c)*I+c*J)/scale; 27 end 28 29 30 f = 1/(1-a-b); 31 C = vech(chol(f*cov(data')))/10; 32 33 if nargout >3 34 par = [C' diag(A)' vech like(A) diag(B)' vech like(B)]; 35 end 36 end APˆ ENDICE C. C ´ ODIGOS MATLAB 57 1function [C, A, B] = par transform(parameters, nob) 2 3dimensoes = nob*(nob+1)/2; 4C = parameters(1:dimensoes); 5 6A = diag(parameters(dimensoes+1: 2*dimensoes)); 7 8A(1,2:3) = parameters(1 + 2*dimensoes:2 + 2*dimensoes); 9A(2,3) = parameters(3 + 2*dimensoes); 10 A(2,1) = A(1,2); 11 A(3,1) = A(1,3); 12 A(3,2) = A(2,3); 13 14 B = diag(parameters(dimensoes *dimensoes+1: (dimensoes+1)*dimensoes)); 15 B(1,2:3) = parameters(1 + (dimensoes+1)*dimensoes:2 + (dimensoes+1)*dimensoes); 16 B(2,3) = parameters(3 + (dimensoes+1)*dimensoes); 17 B(2,1) = B(1,2); 18 B(3,1) = B(1,3); 19 B(3,2) = B(2,3); 20 end APˆ ENDICE C. C ´ ODIGOS MATLAB 64 1%Este script estima e testa a qualidade do ajustamento 2%do modelo usado. 3%Para o seu uso precisa-se carregar os retornos no ambiente de trabalho 4 5if ~ismember({'F','ret'},who) 6clear; 7clc; 8load ret; 9p=1; 10 q=1; 11 [datasym, backcast,T,k] = data transform(ret); 12 13 a = 0.15; 14 b = 0.09; 15 c = 0.07; 16 [~,~,~,x0] = generate par(ret,a,b,c,p,q,k); 17 18 codegen -args {x0,ret,datasym,backcast,p,q,T,k}likelihoods 19 F = @likelihoods mex; 20 end 21 ret(1,:) = ret(1,:) - mean(ret(1,:)); 22 ret(2,:) = ret(2,:) - mean(ret(2,:)); 23 a = 0.03; 24 b = 0.59; 25 c = 0.07; 26 [~,~,~,x0] = generate par(ret,a,b,c,p,q,k); 27 tipo = 'normal'; 28 [parametros, ll, Ht, VCV, scores,hessiano] = mgvech(x0,ret, p, q,tipo,F); 29 30 %Calculo dos erros e criteiros de informacao 31 erros = sqrt(diag(VCV)); 32 erros2 = sqrt(diag(hessiano)); 33 K = length(parametros); 34 [AIC, BIC] = aicbic(ll,K,T); 35 36 37 38 dvec = [0.006101 0.013974 0.02306 0.034965 0.019760 0.051744... 39 0.952451 0.890176 0.911816 ]'; 40 dvec err = [0.016081 0.020129 0.007302 0.006137 0.011796 0.009256 0.00753... 41 0.081961 0.018914]'; 42 43 bekk = [0.057854 0.046324 0.174900 0.153829 -0.053035 0.018800 ... 44 0.254250 0.984434 0.012521 -0.002573 0.940532]'; 45 46 bekk err = [0.008720 0.060460 0.026505 0.015968 0.024387 0.020850 0.028261... 47 0.003495 0.007105 0.008268 0.013902]'; 48 49 nlp ci = [parametros' - 1.96*erros parametros' + 1.96*erros]; 50 nlp ci2 = [parametros' - 1.96*erros2 parametros' + 1.96*erros2]; 51 dvec ci = [dvec - 1.96*dvec err dvec + 1.96*dvec err]; 52 bekk ci = [bekk - 1.96*bekk err bekk + 1.96*bekk err]; 53 54 %% testar a qualidade do ajustamento do modelo 55 load RETORNOS APˆ ENDICE C. C ´ ODIGOS MATLAB 65 56 Hsp500 = zeros(T,1); 57 Hspftse = zeros(T,1); 58 HftseSp500 = zeros(T,1); 59 60 for t=1:T 61 htemp = Ht(:,:,t); 62 Hsp500(t) = htemp(1,1); 63 Hspftse(t) = htemp(2,2); 64 HftseSp500(t) = htemp(1,2); 65 end 66 Z sp500 = (RETORNOS.SP500.ˆ2)./Hsp500; 67 Z ftse = (RETORNOS.FTSE.ˆ2)./Hspftse; 68 69 % Testes sobre a autocorrelacao 70 %corelograma 71 subplot(2,1,1), autocorr(Z ftse) 72 title('Correlograma de residuos padronizados de FTSE') 73 subplot(2,1,2), autocorr(Z sp500) 74 title('Correlograma de residuos padronizados de SP500') 75 76 %testes mais formais: desfazagem 10 15 20 com alpha = 0.05 77 [H sp500,pValue sp500,... 78 Stat sp500,CriticalValue sp500] = lbqtest(Z sp500,[10 15 20]',0.05); 79 [H sp500 pValue sp500 Stat sp500 CriticalValue sp500] 80 81 [H ftse,pValue ftse,... 82 Stat ftse,CriticalValue ftse] = lbqtest(Z ftse,[10 15 20]',0.05); 83 [H ftse pValue ftse Stat ftse CriticalValue ftse] 84 85 %% Testes sobre os efeitos ARCH 86 clc; 87 ZZ sp500 = RETORNOS.SP500./sqrt(Hsp500); 88 ZZ ftse = RETORNOS.FTSE./sqrt(Hspftse); 89 90 [h arch sp500,p arch sp500,... 91 Fstat arch sp500,crit arch sp500] = archtest(ZZ sp500,'lags',2); 92 [h arch sp500,p arch sp500, Fstat arch sp500,crit arch sp500] 93 94 [h arch ftse,p arch ftse,... 95 Fstat arch ftse,crit arch ftse] = archtest(ZZ ftse,'lags',2); 96 [h arch ftse,p arch ftse, Fstat arch ftse,crit arch ftse] 97 98 %% Teste sobre a normalidade dos residuos padronizados 99 clc 100 %qqplots 101 subplot(2,1,1), qqplot(ZZ sp500) 102 subplot(2,1,2), qqplot(ZZ ftse) 103 104 [h sp500 n,p sp500 n,jbstat sp500 n,critval sp500 n] = jbtest(ZZ sp500); 105 [h sp500 n p sp500 n jbstat sp500 n critval sp500 n] 106 107 [h ftse n,p ftse n,jbstat ftse n,critval ftse n] = jbtest(ZZ ftse); 108 [h ftse n p ftse n jbstat ftse n critval ftse n] Apˆendice D C´odigos S-Plus 1#carregar os modulos 2module(finmetrics) 3#Abrir os dados 4 5#formatar o tempo a ser usado 6tempo = timeDate(as.character(anos[,1]), 7+in.format="%d-%m-%y", format="%a %b %d, %Y") 8length(tempo) 9sp500.dados = timeSeries(data = sp500Ftse$Col2,pos = tempo) 10 ftse.dados = timeSeries(data = sp500Ftse$Col3,pos = tempo) 11 12 #formatar alguns elementos dos ativos 13 sp500.dados@title = "Precos de encerramento dos sp500" 14 sp500.dados@units = "Precos em dolar" 15 ftse.dados@title = "Precos de encerramento dos ftse" 16 ftse.dados@units = "Precos em dolar" 17 18 #Calcular os retsp500Ftse 19 sp500.dados.retsp500Ftse = getReturns(sp500.dados,type="continuous", 20 +percentage=T,trim=T) 21 ftse.dados.retsp500Ftse = getReturns(ftse.dados,type="continuous", 22 percentage=T,trim=T) 23 24 #Criar a serie multivariada a partir dos ativos SP500 e FTSE 25 retsp500Ftse = seriesMerge(sp500.dados.retsp500Ftse, 26 +ftse.dados.retsp500Ftse) 27 colIds(retsp500Ftse)<-c("sp500","ftse-20") 28 acf(retsp500Ftseˆ2) 29 30 #-------------------------------------------------------------- 31 # Modelo DVECH 32 #-------------------------------------------------------------- 33 34 #Estimar o modelo DVECH 35 retsp500Ftse.dvec = mgarch(retsp500Ftse~1, ~dvec(1,1), trace=F) 36 retsp500Ftse.dvec 37 38 #Obter as estatisticas do modelo estimado 66 APˆ ENDICE D. C ´ ODIGOS S-PLUS 67 39 summary(retsp500Ftse.dvec) 40 41 #Efectuar alguns teste sobre a qualidade do ajustamento do modelo 42 autocorTest(residuals(retsp500Ftse.dvec, standardize=T)ˆ2,lag=12, bycol=F) 43 archTest(residuals(retsp500Ftse.dvec, standardize=T), lag=12) 44 45 #analise grafica do modelo 46 plot(retsp500Ftse.dvec) 47 48 49 #-------------------------------------------------------------- 50 # Modelo BEKK 51 #-------------------------------------------------------------- 52 53 #Estimar o modelos BEKK 54 retsp500Ftse.bekk = mgarch(retsp500Ftse~1, ~bekk(1,1),trace=F) 55 56 #Obter as estatisticas do modelo estimado 57 summary(retsp500Ftse.bekk) 58 59 #Efectuar alguns teste sobre a qualidade do ajustamento do modelo 60 autocorTest(residuals(retsp500Ftse.bekk, standardize=T)ˆ2,lag=12, bycol=F) 61 archTest(residuals(retsp500Ftse.bekk, standardize=T), lag=12) 62 63 #analise grafica do modelo 64 plot(retsp500Ftse.bekk)