Estimação de parâmetros em modelos estocásticos de estruturas com comportamento dinâmico linear e quasi linear
Full text
Ana Filipa Prior Estima¸c˜ao de parˆametros em modelos estoc´asticos de estruturas com comportamento dinˆamico linear e quasi linear Departamento de Matem´atica da Faculdade de Ciˆencias da Universidade do Porto Janeiro de 2015
Ana Filipa Prior Estima¸c˜ao de parˆametros em modelos estoc´asticos de estruturas com comportamento dinˆamico linear e quasi linear Tese submetida `a Faculdade de Ciˆencias da Universidade do Porto para obten¸c˜ao do grau de Doutor em Matem´atica Aplicada Departamento de Matem´atica da Faculdade de Ciˆencias da Universidade do Porto Janeiro de 2015
Aos meus filhos, por quem tudo fa¸co, por quem tudo sofro, por quem tudo vale a pena. iii
Agradecimentos “Se encontramos um caminho na vida sem obst´aculos provavelmente n˜ao leva a lugar nenhum.” A realiza¸c˜ao de uma tese de Doutoramento ´e um passo importante no percurso de um investigador-professor. Este caminho a ser percorrido ´e longo e ´arduo e seria inexequ´ıvel sem o devido apoio. Aos meus colegas do ISEL, agrade¸co, todo o apoio a n´ıvel profissional, bem como, a for¸ca moral que me fortaleceu e fez permanecer no caminho certo. Ao Pedro Vieira, agrade¸co o seu apoio e o seu trabalho na realiza¸c˜ao dos scripts usados nos estudos com base em simula¸c˜oes nesta tese. ` A Professora Marina Kleptsyna agrade¸co o seu apoio e trabalho colaborativo na Parte I desta tese de Doutoramento. Aos meu Pais, que desde que nasci fazem todos os sacrif´ıcios deste Mundo para que eu seja feliz e concretize todos os meus sonhos, a eles devo tudo, a minha vida, os meus sucessos, a Mulher que sou hoje. Aos meus filhos, lamento pelos momentos em que me ausentei, pela fadiga e temperamento dif´ıcil de alguns dias. Tenho a esperan¸ca que esta tese de Doutoramento lhes sirva de li¸c˜ao de vida, que com grande esfor¸co, dedica¸c˜ao e sacrif´ıcio tudo ´e poss´ıvel, e quero que saibam, que, tal como os meus Pais, eu estarei sempre a seu lado a apoi´a-los para concretizarem os seus sonhos. Ao meu marido agrade¸co o seu apoio, a sua paciˆencia, a sua espera e esperan¸ca num futuro melhor. Aos meu amigos e familiares pe¸co desculpa pela minha ausˆencia e falta de disponibilidade, mas quero que saibam que o seu apoio e incentivo foi importante nas alturas mais dif´ıceis. Quero agradecer em especial `a minha amiga Carla e aos seus Pais, pelo apoio dado a tantos n´ıveis, mas em especial, pelo carinho que sempre me deram e com que sempre me receberam. Ao Centro de Matem´atica da Universidade do Porto agrade¸co todo o apoio facultado na participa¸c˜ao em conferˆencias.
Por ´ultimo, agrade¸co `a pessoa que esteve sempre a meu lado, sem a qual este trabalho n˜ao seria poss´ıvel, com quem desenvolvi uma forte rela¸c˜ao de amizade ao longo deste tempo, que tamb´em fez sacrif´ıcios para me apoiar, que me mostrou outra realidade e modo de pensar e me abriu as portas a outros mundos, a minha orientadora, Professora Paula Milheiro de Oliveira. Este trabalho de investiga¸c˜ao foi feito ao abrigo do programa PROTEC atrav´es do IPL. vi
Resumo Esta disserta¸c˜ao tem como objetivo a investiga¸c˜ao de estimadores de parˆametros de equa¸c˜oes diferenciais estoc´asticas que servem usualmente para modelar o comportamento dinˆamico linear ou quasi linear de estruturas. Para al´em da obten¸c˜ao de estimadores, pretende-se estudar as suas propriedades. A disserta¸c˜ao est´a dividida em 3 partes. Na Parte I abordamos o problema da estima¸c˜ao da matriz de deriva de um modelo linear estoc´astico homog´eneo de dimens˜ao 2ncom coeficientes constantes, observado em tempo cont´ınuo e sendo a matriz de difus˜ao singular. A matriz de deriva ´e uma matriz por blocos figurando nos blocos superiores, n×n, a matriz nula e a matriz identidade e, nos blocos inferiores, n×n, as matrizes M−1KeM−1C, sendo Ma matriz de massa invert´ıvel, Ka matriz de rigidez e Ca matriz de amortecimento. Descrevemos o estimador de m´axima verosimilhan¸ca da matriz de deriva e apresentamos uma demonstra¸c˜ao da sua propriedade da distribui¸c˜ao assint´otica normal com recurso a t´ecnicas da Transformada de Laplace. Estudamos tamb´em a convergˆencia da matriz de covariˆancia deste estimador e explicitamos a matriz de Informa¸c˜ao de Fisher num caso em que se verifica uma condi¸c˜ao de comutatividade da multiplica¸c˜ao de matrizes do modelo. Este estudo ´e acompanhado de simula¸c˜oes que ilustram e complementam os resultados te´oricos obtidos. Na Parte II tratamos o problema de estima¸c˜ao da matriz de deriva de dimens˜ao 2 do mesmo modelo no caso em que ocorre uma mudan¸ca de regime nesta matriz. Esta mudan¸ca de regime acontece quando, num dado instante (conhecido ou desconhecido), o coeficiente de rigidez muda de um determinado valor k1>0 para k2>0 o que nos conduz a um modelo linear por tro¸cos. Considerando as observa¸c˜oes em tempo discreto, descrevemos uma forma de obter estimativas de m´axima verosimilhan¸ca dos parˆametros k1,k2ec(parˆametro de amortecimento) do modelo estoc´astico antes e depois da mudan¸ca de regime. No caso em que a mudan¸ca de regime ocorre num instante desconhecido, este passa a ser outro dos parˆametros a estimar. Mostramos tamb´em como realizar um teste de hip´oteses para a dete¸c˜ao de mudan¸ca de regime, no caso em que se desconhece se efetivamente esta ocorreu. Apresentamos estudos baseados em simula¸c˜oes que ilustram o procedimento e nos permitem analisar a influˆencia dos parˆametros no desempenho dos testes e probabilidades de erro do
xiv
Lista de Tabelas 1.1 Valor m´edio e desvio padr˜ao obtidos a partir das estimativas dos parˆametros kec, para diferentes valores de T(Exemplo 1). . . . . . . . . . . . . . . . . . 36 1.2 Resultados do teste de normalidade de Kolmogorov-Smirnov aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de T, para testar a distribui¸c˜ao assint´otica marginal conforme estabelecida no Teorema 1.4.1 (Exemplo1). .................................... 36 1.3 Resultados do teste de n˜ao correla¸c˜ao para ˆ ke ˆc(Exemplo 1). . . . . . . . . . 36 1.4 Valor m´edio e desvio padr˜ao calculados a partir de estimativas dos parˆametros kec, obtidos para diferentes valores de T(Exemplo 2). . . . . . . . . . . . . 40 1.5 Resultados do teste de normalidade de Kolmogorov-Smirnov aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de T, para testar a distribui¸c˜ao assint´otica estabelecida no Teorema 1.4.1 (Exemplo 2). . . . . . 44 1.6 Valor m´edio das estimativas do parˆametro kobtido para diferentes valores de Te de σ(Exemplo3)................................ 44 1.7 Valor m´edio das estimativas do parˆametro cobtido para diferentes valores de Te de σ(Exemplo3)................................ 45 1.8 Valores m´edios e desvios padr˜ao calculados para as estimativas dos parˆametros das matrizes KeC, obtidos para diferentes valores de T(Exemplo 4). . . . . 49 1.9 Valores m´edios e desvios padr˜ao calculados para as estimativas dos parˆametros das matrizes KeC, obtidas para diferentes valores de T(Exemplo 5). . . . . 51 2.1 Valor m´edio e desvio padr˜ao calculados sobre as estimativas dos parˆametros k,c1ec2, para diferentes valores de T(Exemplo 4). . . . . . . . . . . . . . . 62 xv
3.1 (Exemplo 1) Valor m´edio calculado a partir das estimativas dos parˆametros kt ecobtidos para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer. ....................................... 81 3.2 (Exemplo 1) Valor m´edio calculado a partir das estimativas dos parˆametros kt ecobtidos para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer. ....................................... 85 3.3 M´edias das estimativas dos parˆametros quando o instante de mudan¸ca de regime ´e desconhecido (u= 263 s). ........................ 90 3.4 (Exemplo 3) M´edias das estimativas dos parˆametros para T= 600 squando o instante de mudan¸ca de regime ´e desconhecido (u= 285 s). . . . . . . . . . 96 3.5 (Exemplo 4) M´edias das estimativas dos parˆametros para T= 600 squando o instante de mudan¸ca de regime ´e desconhecido (u= 315 s). . . . . . . . . . 98 3.6 (Exemplo 5) M´edias das estimativas dos parˆametros quando o instante de mudan¸ca de regime ´e desconhecido (u= 329 s). ................. 99 4.1 Valores cr´ıticos para diferentes n´ıveis de significˆancia, α, do teste de hip´oteses (4.4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 4.2 (Exemplo 1) Valor da estat´ıstica de teste obtido sobre uma trajet´oria e valores cr´ıticos para diferentes n´ıveis de significˆancia (o valor simulado foi u= 263 s). 123 4.3 (Exemplo 2) Valor da estat´ıstica de teste obtido sobre uma trajet´oria e valores cr´ıticos para diferentes n´ıveis de significˆancia (o valor simulado do instante de mudan¸ca de regime foi u= 285 s). ........................125 4.4 (Exemplo 3) Valor da estat´ıstica de teste obtido sobre uma trajet´oria e valores cr´ıticos para diferentes n´ıveis de significˆancia (o valor simulado do instante de mudan¸ca de regime foi u= 315 s). ........................127 4.5 (Exemplo 4) Valor da estat´ıstica de teste obtido sobre uma trajet´oria e valores cr´ıticos para diferentes n´ıveis de significˆancia (o valor simulado do instante de mudan¸ca de regime foi u= 329 s). ........................129 5.1 M´edia das estimativas do parˆametro k, obtida para diferentes valores de Te de H(Exemplo1)..................................155 5.2 M´edia das estimativas do parˆametro c, obtida para diferentes valores de Te de H(Exemplo1)..................................156 xvi
5.3 Resultados do teste de normalidade de Shapiro-Wilk aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de TeH= 0.55 (Exemplo 1). ..........................................162 5.4 Resultados do teste de normalidade de Shapiro-Wilk aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de TeH= 0.75 (Exemplo 1). ..........................................162 5.5 Resultados do teste de normalidade de Shapiro-Wilk aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de TeH= 0.95 (Exemplo 1). ..........................................162 xvii
xviii
Lista de Figuras 1 Exemplos de estruturas monitorizadas . . . . . . . . . . . . . . . . . . . . . . 2 1.1 Regi˜oes de confian¸ca para ˆ ke ˆcobtidas para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 1.2 Histogramas para ˆ kobtidos para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva gaussiana aos dados. 38 1.3 Histogramas para ˆcobtidos para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva gaussiana aos dados. . 39 1.4 Regi˜oes de confian¸ca para ˆ ke ˆcobtidas para diferentes valores de T(Exemplo 2): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 1.5 Histogramas para ˆ kobtidos para diferentes valores de T(Exemplo 2): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva Gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 42 xix
1.6 Histogramas para ˆcobtidos para diferentes valores de T(Exemplo 2): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva Gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 43 1.7 Curvas m´edias das estimativas do parˆametro kobtida para diferentes valores de σ(Exemplo3). ................................. 45 1.8 Curvas m´edias das estimativas do parˆametro cobtidas para diferentes valores de σ(Exemplo3). ................................. 46 1.9 RMSE das estimativas do parˆametro kobtidas para diferentes valores de σ (Exemplo3). .................................... 46 1.10 RMSE das estimativas do parˆametro cobtidas para diferentes valores de σ (Exemplo3). .................................... 47 2.1 Histogramas para ˆ kobtidos para diferentes valores de T(Exemplo 4): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. . . . . . . . . . . . . . . . . 63 2.2 Histogramas para ˆc1obtidos para diferentes valores de T(Exemplo 4): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. . . . . . . . . . . . . . . . . 64 2.3 Histogramas para ˆc2obtidos para diferentes valores de T(Exemplo 4): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. . . . . . . . . . . . . . . . . 65 3.1 (Exemplo 1) Histogramas de ˆ ktpara diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer, e curva gaussiana ajustada aos dados: (a)T= 200s, (b)T= 600s. ................................. 82 3.2 (Exemplo 1) Histogramas de ˆcpara diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer, e curva gaussiana ajustada aos dados: (a)T= 200s, (b)T= 600s. ................................. 83 xx
3.3 Regi˜oes de confian¸ca para Ehˆ ktieE[ˆc] para um n´ıvel de confian¸ca de 95% para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer (Exemplo 1): (a)T= 200s, (b)T= 600s. O ponto m´edio da nuvem ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.4 (Exemplo 1) Histogramas de ˆ ktpara diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer, e curva gaussiana ajustada aos dados: (a)T= 200s, (b)T= 600s. ................................. 85 3.5 (Exemplo 1) Histogramas de ˆcpara diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer, e curva gaussiana ajustada aos dados: (a)T= 200s, (b)T= 600s. ................................. 86 3.6 (Exemplo 1) Regi˜oes de confian¸ca para Ehˆ ktieE[ˆc] para um n´ıvel de confian¸ca de 95% para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer: (a)T= 200s, (b)T= 600s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. 87 3.7 (Exemplo 2) Estimativa do coeficiente de rigidez kcalculadas sobre uma trajet´oria ao longo do tempo T(u= 263 s). ................... 91 3.8 (Exemplo 2) Estimativa do coeficiente de amortecimento ccalculadas sobre uma trajet´oria ao longo do tempo T(u= 263 s). ................ 91 3.9 Fun¸c˜ao verosimilhan¸ca para o Exemplo 2 calculada sobre uma das 200 trajet´orias (u= 263 s).................................. 92 3.10 (Exemplo 2) M´edia das estimativas do coeficiente de rigidez calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 263 s). . . . 92 3.11 (Exemplo 2) M´edia das estimativas do coeficiente de amortecimento calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 263 s). 93 3.12 (Exemplo 2) RMSE das estimativas do coeficiente de rigidez (u= 263 s). . . . 94 3.13 (Exemplo 2) RMSE das estimativas do coeficiente de amortecimento (u= 263 s). 94 3.14 (Exemplo 3) M´edia das estimativas do coeficiente de rigidez calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 285 s). . . . 95 3.15 (Exemplo 3) M´edia das estimativas do coeficiente de amortecimento calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 285 s). 96 xxi
3.16 (Exemplo 4) M´edia das estimativas do coeficiente de rigidez calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 315 s). . . . 97 3.17 (Exemplo 4) M´edia das estimativas do coeficiente de amortecimento calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 315 s). 98 3.18 (Exemplo 5) M´edia das estimativas do coeficiente de rigidez calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 329 s). . . . 99 3.19 (Exemplo 5) M´edia das estimativas do coeficiente de amortecimento calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 329 s). 100 4.1 (Exemplo 1) Valores da estat´ıstica de teste para as 200 trajet´orias geradas na simula¸c˜ao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.2 (Exemplo 1) Estimativas do instante de mudan¸ca de regime para as 200 simula¸c˜oes (o valor simulado do instante de mudan¸ca de regime foi u= 263 s). 124 4.3 (Exemplo 2) Valores da estat´ıstica de teste para as 200 trajet´orias geradas na simula¸c˜ao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 4.4 (Exemplo 2) Estimativas do instante de mudan¸ca de regime para as 200 simula¸c˜oes (o valor simulado do instante de mudan¸ca de regime foi u= 285 s). 126 4.5 (Exemplo 3) Valores da estat´ıstica de teste para as 200 trajet´orias geradas na simula¸c˜ao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.6 (Exemplo 3) Estimativas do instante de mudan¸ca de regime para as 200 simula¸c˜oes (o valor simulado do instante de mudan¸ca de regime foi u= 315 s). 128 4.7 (Exemplo 4) Valores da estat´ıstica de teste para as 200 trajet´orias geradas na simula¸c˜ao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 4.8 (Exemplo 4) Estimativas do instante de mudan¸ca de regime para as 200 simula¸c˜oes (o valor simulado do instante de mudan¸ca de regime foi u= 329 s). 130 5.1 M´edia das trajet´orias de ˆ kao longo do tempo, obtida para diferentes valores de H(Exemplo1)..................................156 5.2 M´edia das trajet´orias de ˆcao longo do tempo, obtida para diferentes valores de H(Exemplo1)..................................157 5.3 RMSE de ˆ kao longo do tempo, para diferentes valores de H(Exemplo 1). . . 157 5.4 RMSE de ˆcao longo do tempo, para diferentes valores de H(Exemplo 1). . . 158 xxii
xxiii 5.5 Nuvens de pontos e regi˜oes de confian¸ca para ˆ ke ˆccom H= 0.55,obtidas para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. . . . . . . . . . . . . . . . 159 5.6 Nuvens de pontos e regi˜oes de confian¸ca para ˆ ke ˆccom H= 0.75,obtidas para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. . . . . . . . . . . . . . . . 160 5.7 Nuvens de pontos e regi˜oes de confian¸ca para ˆ ke ˆccom H= 0.95,obtidas para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. . . . . . . . . . . . . . . . 161 5.8 Histogramas para ˆ kcom H= 0.55 obtidos para diferentes valores de T (Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 163 5.9 Histogramas para ˆ kcom H= 0.75 obtidos para diferentes valores de T (Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 164 5.10 Histogramas para ˆ kcom H= 0.95 obtidos para diferentes valores de T (Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 165 5.11 Histogramas para ˆccom H= 0.55 obtidos para diferentes valores de T (Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 166 5.12 Histogramas para ˆccom H= 0.75 obtidos para diferentes valores de T (Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 167 5.13 Histogramas para ˆccom H= 0.95 obtidos para diferentes valores de T (Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a cheio representa o ajustamento de uma curva Gaussiana aos dados. 168
4CAP´ ITULO 0. INTRODUC¸ ˜ AO vista das aplica¸c˜oes, o modelo em dimens˜ao 2napresenta um maior interesse pr´atico, dada a complexidade de estruturas que representa, contendo nn´os. A utiliza¸c˜ao do modelo (2) num outro ˆambito de aplica¸c˜ao, o dos circuitos el´etricos, pode ser consultado em [2]. De facto, o modelo que surge em alguns problemas de circuitos el´etricos ´e matematicamente semelhante ao da mecˆanica estrutural, embora tendo obviamente diferente interpreta¸c˜ao. Nesta disserta¸c˜ao, por quest˜oes de clareza e simplicidade opt´amos frequentemente por usar a terminologia da mecˆanica estrutural (matrizes de massa, matrizes de rigidez e matrizes de amortecimento) que nos ´e mais familiar do que a terminologia dos circuitos el´etricos. O processo estoc´astico que resulta do modelo dinˆamico estoc´astico linear com coeficientes constantes (modelo (2)) ´e um processo de difus˜ao habitualmente designado na literatura por processo de Ornstein-Uhlenbeck (cf [2], por exemplo) e foi objeto de investiga¸c˜ao por diversos autores, devido essencialmente `a sua capacidade para modelar fen´omenos f´ısicos, apesar da sua simplicidade. Trabalhos como [3], [4] e [5] analisam o problema de estima¸c˜ao de parˆametros do processo de Ornstein-Uhlenbeck atrav´es da constru¸c˜ao da fun¸c˜ao de verosimilhan¸ca. Apesar destes autores se situarem em cen´arios que apresentam importantes restri¸c˜oes, revelam-se particularmente ´uteis pelas metodologias que desenvolveram. De facto, o m´etodo de m´axima verosimilhan¸ca tem sido sugerido na literatura para resolver este tipo de problemas de estima¸c˜ao de parˆametros. Neste sentido, a revis˜ao da metodologia de determina¸c˜ao da derivada de Radon-Nikodym para a constru¸c˜ao da fun¸c˜ao de verosimilhan¸ca assume particular importˆancia na obten¸c˜ao dos estimadores. Esta t´ecnica de constru¸c˜ao da fun¸c˜ao de verosimilhan¸ca est´a documentada por exemplo em [2], [10] e [67]. No entanto somente em [2] e [67] ´e contemplado o caso da derivada de Radon-Nikodym num contexto mais geral duma matriz de deriva degenerada, que corresponde ao caso que nos interessa diretamente nesta tese (modelo (2)). Esta t´ecnica permite, utilizando pseudo-inversas da matriz B, construir os estimadores em contextos muito gerais. No entanto, no que se refere `a dedu¸c˜ao das propriedades dos estimadores os resultados te´oricos apresentados em [2] restringem-se ao caso n˜ao degenerado. Ainda no que diz respeito a esta metodologia s˜ao igualmente de salientar os trabalhos [25], [61], [63], [82] (onde a matriz B´e considerada uma matriz definida positiva) e o trabalho [47] (onde a matriz de difus˜ao ´e singular mas o modelo em estudo n˜ao ´e o modelo estoc´astico linear com coeficientes constantes que nos propomos
5 estudar). Grande parte da literatura dedicada `a estima¸c˜ao de parˆametros de processos de difus˜ao aborda apenas as difus˜oes unidimensionais (ver, por exemplo, [4] ou [19]). A investiga¸c˜ao no ˆambito da estima¸c˜ao de parˆametros de processos de difus˜ao no caso multidimensional e das propriedades dos estimadores ´e bastante mais recente. Neste contexto, podemos salientar os resultados obtidos em [2], [9], [50] e [66] respeitantes `a estima¸c˜ao de parˆametros de processo de difus˜ao no caso multidimensional e `as propriedades dos estimadores obtidos. ´ E sabido que, no estudo do problema da estima¸c˜ao dos parˆametros A e B do modelo (2), considerando observa¸c˜oes em tempo cont´ınuo, a estima¸c˜ao da matriz Brecorrendo `a F´ormula de Itˆo n˜ao oferece dificuldade. Por conseguinte, apenas a estima¸c˜ao da matriz A constitui um problema que merece investiga¸c˜ao. O estudo da estima¸c˜ao da matriz Bapenas se torna interessante no caso em que as observa¸c˜oes ocorrem em tempo discreto (veja-se por exemplo [29]), o que n˜ao corresponde ao nosso ponto de partida nesta tese. Outros trabalhos publicados sobre a estima¸c˜ao de parˆametros em modelos em tempo cont´ınuo traduzidos por equa¸c˜oes diferenciais estoc´asticas, utilizam abordagens distintas, optando por discretiz´a-los, aproximando-os por processos em tempo discreto, ou seja, por s´eries temporais ([8]). Essa discretiza¸c˜ao pode ser feita pelos esquemas habituais, estoc´asticos, existentes na literatura (ver por exemplo [59]). Uma vez investigado o problema de estima¸c˜ao nos modelos lineares estamos em condi¸c˜oes de fazer uma incurs˜ao pelos modelos quasi lineares, mais propriamente tratando modelos lineares por tro¸cos. Vamos investigar o problema de estima¸c˜ao da matriz de deriva, A, do modelo (1), sujeito a observa¸c˜oes em tempo discreto quando, ap´os um certo per´ıodo de tempo, ocorre uma mudan¸ca de regime nesta matriz e por consequˆencia na dinˆamica do processo. Referimos que, em engenharia, reveste-se de particular importˆancia a identifica¸c˜ao de altera¸c˜oes no parˆametro de rigidez da matriz de deriva, em particular uma diminui¸c˜ao do parˆametro de rigidez, o que representa a existˆencia de um dano estrutural. Encontramos exemplos de modelos com mudan¸ca de regime noutras ´areas como a matem´atica financeira ou a biologia, considerando a possibilidade de ocorrerem v´arias mudan¸cas de regime no intervalo de tempo considerado (ver [31] e [48]). No entanto, no caso da engenharia, esta possibilidade de haver v´arias mudan¸ca de regime n˜ao faz muito sentido, pois, uma vez ocorrida a mudan¸ca de regime, a estrutura encontra-se danificada. Mais uma vez
6CAP´ ITULO 0. INTRODUC¸ ˜ AO o m´etodo da m´axima verosimilhan¸ca ´e usado para estimar os parˆametros (de rigidez e de amortecimento) que constituem a matriz de deriva A, antes e depois da mudan¸ca de regime ocorrer. Usamos as probabilidades de transi¸c˜ao condicionadas dos estados para construir a fun¸c˜ao de verosimilhan¸ca. No entanto, devido `a complexidade da fun¸c˜ao de verosimilhan¸ca o estimador de m´axima verosimilhan¸ca dos parˆametros ´e obtido apenas por otimiza¸c˜ao num´erica desta fun¸c˜ao. Por este motivo, opt´amos por abordar este problema considerando, numa primeira etapa, as observa¸c˜oes do processo em tempo cont´ınuo. A fun¸c˜ao de verosimilhan¸ca ´e constru´ıda usando a derivada de Radon-Nikodym como no caso do modelo linear, com as necess´arias adapta¸c˜oes. Esta abordagem permite explicitar o estimador de m´axima verosimilhan¸ca dos parˆametros e obtem-se uma aproxima¸c˜ao ao estimador de m´axima verosimilhan¸ca em tempo discreto discretizando este estimador. Uma abordagem semelhante j´a foi antes utilizada noutro contexto (ver [21] e [93]). Devido `as vantagens com que nos depar´amos ao comparar as duas abordagens tamb´em utiliz´amos esta ´ultima abordagem quando consider´amos que o instante em que ocorre a mudan¸ca de regime ´e desconhecido e mostramos como aplicar o m´etodo da m´axima verosimilhan¸ca neste caso, pois o instante de mudan¸ca de regime passa a ser um dos parˆametros a estimar. Em modelos em que h´a suspeita de ocorrˆencia de mudan¸ca de regime, o estudo da estima¸c˜ao dos parˆametros destes modelos deve ser precedido da realiza¸c˜ao de um teste de hip´oteses para dete¸c˜ao da mudan¸ca de regime e, por este motivo, o problema de realiza¸c˜ao do referido teste de hip´oteses ´e abordado nesta disserta¸c˜ao. Podemos ver em [11], [28] e [61] uma revis˜ao e actualiza¸c˜ao dos resultados existentes relativos a testes de hip´oteses do tipo ”change point”. No entanto, estamos particularmente interessados na t´ecnica usada em [32] para a constru¸c˜ao de um teste de hip´oteses para a dete¸c˜ao de mudan¸ca de regime em modelos autoregressivos em tempo discreto e na adapta¸c˜ao que foi feita desta t´ecnica em [34] a modelos observados em tempo cont´ınuo mas com EDEs revers´ıveis para uma m´edia peri´odica. O conte´udo dessas obras inspira a abordagem que adot´amos no nosso problema em concreto, pois, conforme referimos, uma abordagem em tempo cont´ınuo envolvendo uma discretiza¸c˜ao mais tardia do procedimento ´e pass´ıvel de fornecer bons resultados, permitindo uma not´oria redu¸c˜ao do esfor¸co computacional. Sendo assim, vamos seguir a mesma ideia na constru¸c˜ao de um teste de hip´oteses para a ocorrˆencia de mudan¸ca de regime, mantendo o problema inicialmente como se se baseasse em observa¸c˜oes em tempo cont´ınuo para obter a estat´ıstica de teste e s´o depois entramos na fase de discretiza¸c˜ao. A estat´ıstica de teste
7 ´e obtida atrav´es da raz˜ao das verosimilhan¸cas e prova-se que o supremo da estat´ıstica de teste converge em distribui¸c˜ao para uma ponte browniana bidimensional. Alguns autores, como por exemplo [34], s˜ao de opini˜ao que a obten¸c˜ao da convergˆencia da estat´ıstica de teste para a distribui¸c˜ao de uma ponte browniana bidimensional ´e insuficiente para a sua utiliza¸c˜ao em termos pr´aticos, porque consideram ser necess´ario, nesse caso, simular a referida ponte browniana para determinar os valores cr´ıticos do teste. No entanto, n˜ao partilhamos dessa opini˜ao. Encontramos diversos autores (ver [40], [45], [46] e [51], por exemplo) que se dedicaram a estudar m´etodos num´ericos para o c´alculo da distribui¸c˜ao das pontes brownianas e construiram tabelas para a referida distribui¸c˜ao suficientemente precisas para nos permitirem aplicar o teste. Apresentamos um exemplo de aplica¸c˜ao do teste de hip´oteses referido acima e ilustramos o procedimento a ser seguido num estudo com base em simula¸c˜oes. Por ´ultimo tratamos o problema de estima¸c˜ao da matriz de deriva, A, do modelo (1) com base em observa¸c˜oes de {Xt}em tempo cont´ınuo, num intervalo [0, T], agora com o processo de Wiener {Wt}substitu´ıdo por um movimento browniano fracion´ario (MBF) de dimens˜ao 2, WH t= WH 1,t WH 2,t , onde {WH 1,t}e{WH 2,t}s˜ao movimentos brownianos fracion´arios indepedentes com o mesmo parˆametro de Hurst, H, com H∈(1/2,1). Este modelo que consideramos estende o modelo (1) bidimensional. Modelos em que intervˆem ru´ıdos fracion´arios ainda n˜ao s˜ao muito usados em engenharia. No entanto, muitos problemas, em engenharia, fazem intervir processos com mem´oria, e, por este motivo, o estudo do problema de estima¸c˜ao dos parˆametros de modelos lineares perturbados por movimentos brownianos fracion´arios (MBF) tem potencial de aplica¸c˜ao em diversas ´areas da engenharia (ver [14], [39] e [70]). Para este novo modelo n˜ao ´e poss´ıvel aplicar diretamente as t´ecnicas de inferˆencia estat´ıstica a que recorremos no modelo (1), visto que movimentos brownianos fracion´arios n˜ao s˜ao martingalas. Por este mesmo motivo, a teoria cl´assica da integra¸c˜ao estoc´astica n˜ao pode ser aplicada, pelo menos diretamente, sendo esta a principal dificuldade em lidar com este modelo. No entanto, vamos conseguir aplicar o m´etodo da m´axima verosimilhan¸ca para obter o EMV do parˆametro de rigidez e do parˆametro de amortecimento que constituem a matriz de deriva, A. De facto, vamos mostrar que uma transforma¸c˜ao adequada deste modelo, recorrendo `a chamada martingala fundamental introduzida em [79], permite usar o Teorema de Girsanov usual para construir a derivada de Radon-Nikodym.
8CAP´ ITULO 0. INTRODUC¸ ˜ AO Por consequˆencia conseguimos obter a express˜ao da fun¸c˜ao de verosimilhan¸ca associada ao problema de estima¸c˜ao de parˆametros. Vamos seguir o percurso sugerido nos trabalhos [54], [58] e [79] na resolu¸c˜ao do problema de estima¸c˜ao dos parˆametros do modelo. No trabalho que seguidamente se apresenta decidimos, por op¸c˜ao, n˜ao apresentar um cap´ıtulo de resultados te´oricos gerais, j´a conhecidos e necess´arios para o trabalho desenvolvido. Limitar-nos-emos a relembrar que, para ser poss´ıvel atingir os objetivos centrais desta disserta¸c˜ao, s˜ao particularmente importantes a teoria das equa¸c˜oes diferenciais estoc´asticas, as propriedades dos processos estoc´asticos envolvidos nos modelos em estudo (vejam-se os resultados apresentados em [7], [49], [86]), bem como a teoria da inferˆencia estat´ıstica para processos estoc´asticos aplicada na pesquisa dos estimadores dos parˆametros das equa¸c˜oes diferenciais estoc´asticas ([10], [61], [67], [82]). Nas teorias referidas anteriormente revelam-se fundamentais alguns resultados da teoria de martingalas a utilizar no c´alculo estoc´astico e alguns procedimentos gen´ericos e t´ecnicas de estima¸c˜ao de parˆametros de equa¸c˜oes diferenciais estoc´asticas. No que respeita `a estacionaridade e ergodicidade dos processos de difus˜ao, s˜ao particularmente ´uteis os resultados apresentados em [7], [49] e [82]. Nesta linha de pensamento, esta disserta¸c˜ao est´a dividida em trˆes partes e um cap´ıtulo final de conclus˜oes. A Parte I ´e composta por dois cap´ıtulos (Cap´ıtulos 1 e 2) e ´e dedicada ao problema de estima¸c˜ao da matriz de deriva Ado modelo (2) considerando as observa¸c˜oes em tempo cont´ınuo. Obtemos o EMV da matriz de deriva Ae conseguimos demonstrar a propriedade LAN do EMV (Cap´ıtulo 1), o que representa um contributo original, uma vez que a literatura existente impunha restri¸c˜oes `a matriz de difus˜ao, B, mas tamb´em pela t´ecnica utilizada na sua demonstra¸c˜ao com recurso a Transformadas de Laplace. O resultado te´orico obtido ´e refor¸cado com o estudo com base em simula¸c˜oes que ´e realizado. Exploramos em que casos ´e poss´ıvel explicitar a matriz de Informa¸c˜ao de Fisher e mostramos que, se os blocos inferiores da matriz de deriva Ado modelo (2) (isto ´e, as matrizes M−1KeM−1C, onde M, KeCrepresentando, respetivamente, as matrizes de massa, rigidez e amortecimento da EDE que traduz o modelo mecˆanico) e a matriz Σ verificarem a propriedade de comutatividade da multiplica¸c˜ao, ent˜ao uma express˜ao expl´ıcita da matriz de covariˆancia do EMV pode ser deduzida facilmente. Tamb´em abordamos a quest˜ao da estima¸c˜ao da matriz de deriva, A, do modelo (2), em casos em que este modelo apresenta certas particularidades (Cap´ıtulo 2), como por exemplo, as matrizes AeB1 2serem sim´etricas e/ou diagonais, particularidades
9 estas que levam a diferen¸cas not´orias no problema de estima¸c˜ao. No caso em que as matrizes AeB1 2do modelo (2) s˜ao diagonais conseguimos demonstrar a propriedade LAN do EMV. No caso em que a matriz K´e diagonal e a matriz C´e sim´etrica, no modelo (2) em dimens˜ao 4, obtemos o EMV de KeCe realizamos um estudo com base em simula¸c˜oes numa tentativa de analisar o seu comportamento assint´otico, mas n˜ao estabelecemos teoricamente a propriedade LAN do estimador. A Parte II ´e composta por dois cap´ıtulos (Cap´ıtulos 3 e 4) e ´e dedicada ao problema de estima¸c˜ao da matriz de deriva, A, do modelo (1), sujeito a observa¸c˜oes em tempo discreto, quando ocorre uma mudan¸ca de regime nesta matriz e, por consequˆencia, na dinˆamica do processo. Assumimos que, num dado instante (que pode ser conhecido ou desconhecido), o modelo passa de um regime para o outro. Esta mudan¸ca ocorre na matriz de deriva que caracteriza o modelo, conforme se especifica a seguir: o modelo ´e representado pela equa¸c˜ao diferencial estoc´astica (1), onde {Wt}´e um processo de Wiener bidimensional e AteB1 2s˜ao matrizes quadradas 2 ×2; a matriz de deriva Ate a matriz de difus˜ao Bs˜ao as matrizes At= 0 1 −kt m−c m eB= 0 0 0σ2 ,(3) onde kt´e o coeficiente de rigidez, c > 0 ´e o coeficiente de amortecimento do modelo e σ > 0 ´e o desvio-padr˜ao da pertuba¸c˜ao; a mudan¸ca de regime ocorre quando, num dado instante (que pode ser conhecido ou desconhecido), o coeficiente de rigidez ktpassa de um determinado valor k1para outro valor k2, sendo k1, k2>0. A EDE (1) – (3) modela o movimento vibrat´orio de estruturas sujeitas a a¸c˜oes aleat´orias em que h´a mudan¸ca de regime do comportamento dinˆamico da estrutura. Em particular, descreve o cen´ario de uma estrutura que sofreu uma diminui¸c˜ao de rigidez ap´os um determinado per´ıodo de tempo de utiliza¸c˜ao. Considerando as observa¸c˜oes em tempo discreto, descrevemos como obter estimativas de m´axima verosimilhan¸ca dos parˆametros k1,k2ecdo modelo estoc´astico, antes e depois da mudan¸ca de regime. No caso em que a mudan¸ca de regime ocorre num instante desconhecido, este tamb´em passa a ser um dos parˆametros a estimar. Neste caso, para al´em de estabelecermos procedimentos para obter as estimativas de m´axima verosimilhan¸ca dos parˆametros k1,k2,c, obtemos ainda as estimativas do instante de mudan¸ca de regime. Para as situa¸c˜oes em que se desconhece se efetivamente ocorreu uma mudan¸ca de regime no modelo em estudo, construimos um teste de hip´oteses para a dete¸c˜ao de mudan¸ca de regime.
10 CAP´ ITULO 0. INTRODUC¸ ˜ AO Para isso obtemos uma representa¸c˜ao expl´ıcita da estat´ıstica de teste e a sua distribui¸c˜ao assint´otica. Apresentamos, na disserta¸c˜ao, um estudo com base em simula¸c˜oes para ilustrar o procedimento de teste e analisamos os resultados obtidos. O nosso contributo mais importante nesta Parte II reside na abordagem original realizada ao problema de estima¸c˜ao que permitiu obter o EMV dos parˆametros desconhecidos e uma sua aproxima¸c˜ao com boa precis˜ao e custos computacionais n˜ao muito elevados, bem como, na abordagem seguida na constru¸c˜ao do teste de hip´oteses, que permitiu adaptar um teste algo semelhante existente na literatura para modelos em dimens˜ao 1 (ver [34]), sendo de salientar que o modelo em estudo nesta tese ´e um modelo estoc´astico de dimens˜ao 2 com 4 parˆametros desconhecidos. A Parte III ´e composta por apenas um cap´ıtulo (Cap´ıtulo 5) e ´e dedicada ao problema da estima¸c˜ao da matriz A, matriz de deriva, do modelo (1) considerando as observa¸c˜oes em tempo cont´ınuo, mas agora com o processo de Wiener {Wt}substitu´ıdo por um movimento browniano fracion´ario (MBF) de dimens˜ao 2 com componentes independentes. O modelo que consideramos nesta Parte III estende o modelo de dimens˜ao 2 apresentado na Parte I desta disserta¸c˜ao. A forma que adot´amos para tratar o problema de estima¸c˜ao em modelos lineares perturbados por MBF passa por realizar uma transforma¸c˜ao do movimento browniano fracion´ario, visto que movimentos brownianos fracion´arios n˜ao s˜ao martingalas. A referida transforma¸c˜ao permite obter um modelo equivalente em termos estat´ısticos sobre o qual ´e aplic´avel o c´alculo estoc´astico de Itˆo. Conseguimos explicitar a nova fun¸c˜ao de logverosimilhan¸ca e obter o EMV depois de descrito o modelo transformado. Para al´em disso, mostramos a centricidade do estimador. Tamb´em obtemos uma representa¸c˜ao da matriz de covariˆancia do processo transformado que mostramos depender de fun¸c˜oes de Bessel matriciais do tipo 1, suscet´ıvel de ser usada no futuro na determina¸c˜ao da Transformada de Laplace constru´ıda com base nas observa¸c˜oes do processo para, `a semelhan¸ca do Cap´ıtulo 1, conseguir estabelecer a propriedade LAN. Implement´amos o EMV em R e realiz´amos um estudo com base em simula¸c˜oes com o objetivo de analisar o seu comportamento em fun¸c˜ao do tempo e do parˆametro de Hurst. Constat´amos fazer sentido que o EMV verifique a propriedade LAN no exemplo de aplica¸c˜ao estudado e, portanto, fazer sentido a conjetura da propriedade LAN ser extensiva aos modelos lineares da forma (2) com ru´ıdos modelados por MBF. Este Cap´ıtulo 5 ´e um contributo original para a investiga¸c˜ao do problema da estima¸c˜ao de parˆametros em modelos lineares estoc´asticos multidimensionais com movimento browniano fracion´ario, observados em tempo cont´ınuo, e representa o ponto de partida do
11 trabalho futuro resultante desta tese. ´ E importante salientar que o trabalho computacional de programa¸c˜ao e de simula¸c˜oes subjacente `as investiga¸c˜oes desenvolvidas foi realizado pelos intervenientes nesta disserta¸c˜ao. O software utilizado corresponde ao programa R em simultˆaneo com o projeto ”yuima”desenvolvido para inferˆencia e simula¸c˜ao de equa¸c˜oes diferenciais estoc´asticas. A disserta¸c˜ao termina com um cap´ıtulo (Cap´ıtulo 6) contendo as principais conclus˜oes e os coment´arios finais sugeridos pelo trabalho que realiz´amos indicando novas perspetivas de an´alise.
12 CAP´ ITULO 0. INTRODUC¸ ˜ AO
Parte I Modelos estoc´asticos com comportamento dinˆamico linear 13
20 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO ´e poss´ıvel explicitar a matriz de Informa¸c˜ao de Fisher. Na ´ultima sec¸c˜ao deste Cap´ıtulo (Sec¸c˜ao 1.6) realizamos um estudo com base em simula¸c˜oes para refor¸car os resultados te´oricos obtidos e colher alguns detalhes. 1.2 Descri¸c˜ao do problema Seja (Ω,F,{Ft}t, P), um espa¸co de probabilidade filtrado. Consideremos a EDE dXt=AXtdt +B1 2dWt,(1.3) onde {Xt}representa um processo em R2ne{Wt}´e um processo de Wiener com valores em R2n. As matrizes AeB1 2s˜ao definidas por: A= 0Id −M−1K−M−1C eB1 2= 0 0 0 Σ1 2 ,(1.4) sendo Σ uma matriz real definida positiva de dimens˜ao n×n. As matrizes M,KeCs˜ao matrizes reais de dimens˜ao n×n, sendo a matriz Minvert´ıvel. Supomos que os valores pr´oprios da matriz As˜ao complexos com parte real negativa sendo os processos {Wi,t},com i= 1, ..., 2nprocessos de Wiener independentes no espa¸co de probabilidade (Ω,F, P). Consideremos ainda a condi¸c˜ao inicial X0=x0, com x0um vetor conhecido em R2n. Assumimos que a matriz M´e conhecida e as matrizes KeCdesconhecidas. Admite-se que o processo {Xt}´e observado no intervalo [0, T ],sem perda de generalidade. Consideramos portanto o problema de estimar as matrizes KeCdo modelo (1.3) – (1.4), o mesmo ´e dizer que o nosso problema consiste em estimar a matriz de deriva, A, no instante T, o que ´e equivalente a estimar o vetor coluna θ= (θi)∈Rn(2n)×1, que representa a concatena¸c˜ao das nlinhas inferiores da matriz A, linha por linha, com base nas observa¸c˜oes do processo {Xt}, em tempo cont´ınuo, no intervalo [0, T]. No que respeita `a estima¸c˜ao da matriz Σ, como referido na introdu¸c˜ao, em tempo cont´ınuo, uma estimativa pode ser calculada em qualquer intervalo de tempo finito com probabilidade um. Por este motivo, sem perda de generalidade, consideramos a matriz Σ conhecida, nesta tese.
1.3. ESTIMADOR DE M ´ AXIMA VEROSIMILHANC¸A 21 Note-se que, de facto, podemos escolher como representa¸c˜ao alternativa do processo {Xt}a equa¸c˜ao dXt=AXtdt + 0 Σ1 2 dWt, em que {Wt}´e um processo de Wiener com valores em Rn. No entanto, preferimos a representa¸c˜ao (1.3) – (1.4), de modo a evidenciar a singularidade da matriz Be da´ı a exclus˜ao deste modelo da teoria geral. 1.3 Estimador de m´axima verosimilhan¸ca O EMV de Aou de θ, como preferirmos, obt´em-se, como ´e natural, maximizando a fun¸c˜ao de log-verosimilhan¸ca, que ´e, no caso das EDEs, constru´ıda `a custa da derivada de RadonNikodym da medida gerada pelo processo {Xt}com respeito `a medida de Wiener (cf. [10] referˆencia a esta teoria geral). Por uma quest˜ao de simplicidade e de clareza na exposi¸c˜ao vamos come¸car por apresentar o EMV para o modelo de dimens˜ao 2, o que corresponde a ter o caso n= 1. 1.3.1 O caso n= 1 Neste caso, o modelo (1.3) tem como matriz de deriva e matriz de difus˜ao, respetivamente, as seguintes matrizes: A= 0 1 −k m−c m (1.5) e B1 2= 0 0 0σ .(1.6) Na terminologia da mecˆanica das vibra¸c˜oes e da dinˆamica estrutural, mrepresenta a massa (m > 0), ko coeficiente de rigidez (k > 0) , c o coeficiente de amortecimento (c > 0) e σo desvio padr˜ao da perturba¸c˜ao (σ > 0). A fun¸c˜ao de verosimilhan¸ca ´e constru´ıda usando a derivada de Radon-Nikodym dPX dPW(ver [10]), sendo o EMV de θobtido por maximiza¸c˜ao da fun¸c˜ao de log-verosimilhan¸ca. Por outras palavras, a derivada de Radon-Nikodym ´e usada para obter os estimadores de m´axima
22 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO verosimilhan¸ca dos coeficientes de rigidez e amortecimento da equa¸c˜ao diferencial estoc´astica (1.3) em dimens˜ao 2, com AeB1 2definidas em (1.5) e (1.6). Ora, essa derivada faz intervir a inversa da matriz B=B1 2B1 2∗= 0 0 0σ2 .Devido `a particularidade de Bser uma matriz degenerada, ser´a usada uma pseudo-inversa de B, mais especificamente, a inversa generalizada de Moore-Penrose que tem a importante propriedade de ser ´unica. O mesmo ´e dizer que, a pseudo-inversa da matriz B´e a matriz B+= 0 0 01 σ2 (ver [89], Teorema 5.6). A derivada de Radon-Nikodym da medida gerada pelo processo {Xt}, solu¸c˜ao da equa¸c˜ao diferencial estoc´astica (1.3) em dimens˜ao 2 com AeB1 2definidas por (1.5) e (1.6) com respeito `a medida de Wiener ´e dada por (ver [2], Teorema 2): dPX dPW (Xt) = exp tr B+A T Z 0 XtdX∗ t−1 2A∗B+A T Z 0 XtX∗ tdt .(1.7) A proposi¸c˜ao seguinte estabelece a express˜ao do EMV, que daqui se deduz. Proposi¸c˜ao 1.3.1 O estimador de m´axima verosimilhan¸ca para θ= −k m −c m ´e ˆ θT= − ∧ k m −∧ c m = T Z 0 XtX∗ tdt −1T Z 0 XtdX2,t . Demonstra¸c˜ao Atendendo `a express˜ao (1.7), a fun¸c˜ao de log-verosimilhan¸ca L(k, c) = ln dPX dPW (Xt) ´e dada por: L(k, c) = −k mσ2 T Z 0 X1,tdX2,t −c mσ2 T Z 0 X2,tdX2,t −1 2m2σ2 k2 T Z 0 X2 1,tdt + 2ck T Z 0 X1,tX2,tdt +c2 T Z 0 X2 2,tdt . Determinar o estimador de m´axima verosimilhan¸ca de θequivale, como ´e bem conhecido, a resolver o sistema: ∂L ∂k = 0 ∂L ∂c = 0 .
1.3. ESTIMADOR DE M ´ AXIMA VEROSIMILHANC¸A 23 Neste caso, trata-se do sistema: −1 mσ2 T Z 0 X1,tdX2,t −k m2σ2 T Z 0 X2 1,tdt −c m2σ2 T Z 0 X1,tX2,tdt = 0 −1 mσ2 T Z 0 X2,tdX2,t −k m2σ2 T Z 0 X1,tX2,tdt −c m2σ2 T Z 0 X2 2,tdt = 0 ⇐⇒ m T Z 0 X1,tdX2,t +k T Z 0 X2 1,tdt +c T Z 0 X1,tX2,tdt = 0 m T Z 0 X2,tdX2,t +k T Z 0 X1,tX2,tdt +c T Z 0 X2 2,tdt = 0 . Na forma matricial, podemos escrever T Z 0 X2 1,tdt T Z 0 X1,tX2,tdt T Z 0 X1,tX2,tdt T Z 0 X2 2,tdt k m c m + T Z 0 X1,tdX2,t T Z 0 X2,tdX2,t = 0 0 . Obt´em-se o estimador de m´axima verosimilhan¸ca: − ∧ k m −∧ c m = T Z 0 X2 1,tdt T Z 0 X1,tX2,tdt T Z 0 X1,tX2,tdt T Z 0 X2 2,tdt −1 T Z 0 X1,tdX2,t T Z 0 X2,tdX2,t = T Z 0 XtX∗ tdt −1T Z 0 XtdX2,t o que conclui a demonstra¸c˜ao. Ap´os esta breve exposi¸c˜ao respeitante ao modelo (1.3) – (1.4) em dimens˜ao 2, poderemos mais facilmente compreender o estudo deste modelo em dimens˜ao 2n. Come¸caremos pela dedu¸c˜ao do EMV de θe o principal objetivo ser´a a demonstra¸c˜ao da propriedade LAN para o EMV. 1.3.2 O caso n≥2 Como referido, o EMV associado ao modelo (1.3) – (1.4) pode ser obtido, como usualmente, por maximiza¸c˜ao da fun¸c˜ao de log-verosimilhan¸ca, a que chamaremos L(θ) (ver [10]), onde
24 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO L(θ) = ln dPX dPW (Xt). Tem-se ent˜ao: L(θ) = tr B+A T Z 0 XtdX∗ t−1 2A∗B+A T Z 0 XtX∗ tdt . Repetindo um racioc´ınio an´alogo ao realizado na Proposi¸c˜ao 1.3.1, obtemos o EMV associado ao modelo (1.3) – (1.4) resolvendo o sistema: ∂L ∂k1 = 0 . . . ∂L ∂kn = 0 ∂L ∂c1 = 0 . . . ∂L ∂cn = 0 . Como alternativa, para evitar todos os c´alculos que da´ı advˆem, podemos aplicar o resultado obtido em [60] ao modelo (1.3) – (1.4) e obtemos o EMV de Adado por: ˆ A= 0Id ˆ An:2n , com ˆ An:2n=ZT 0 Xtd(Xn+1,t Xn+2,t ···X2n,t)∗·ZT 0 XtX∗ tdt−1 , onde An:2nrepresenta a submatriz das n´ultimas linhas de Acontendo os parˆametros desconhecidos. Fazendo uso da ´algebra de matrizes, verificamos que o EMV de θ´e dado
1.3. ESTIMADOR DE M ´ AXIMA VEROSIMILHANC¸A 25 por: ˆ θT=Diag ZT 0 XtX∗ tdt−1 n blocos! ZT 0 XtdXn+1,t ZT 0 XtdXn+2,t . . . ZT 0 XtdX2n,t = Idn×n⊗ZT 0 XtX∗ tdt−1!· ZT 0 XtdXn+1,t ZT 0 XtdXn+2,t . . . ZT 0 XtdX2n,t (1.8) Centricidade, consistˆencia e eficiˆencia de ˆ θT A partir de (1.3) podemos escrever que ZT 0 Xtd[Xn+1,t Xn+2,t ···X2n,t] = ZT 0 XtX∗ tdtA∗ n:2n +ZT 0 Xtd[Wn+1,t Wn+2,t ···W2n,t]Σ1 2 ∗ , isto ´e ZT 0 XtdXn+1,t ZT 0 XtdXn+2,t . . . ZT 0 XtdX2n,t =Idn×n⊗ZT 0 XtX∗ tdtθ+Σ1 2⊗Id2n×2n· ZT 0 XtdWn+1,t ZT 0 XtdWn+2,t . . . ZT 0 XtdW2n,t . Assim, usando a express˜ao (1.8), depois de alguma manipula¸c˜ao alg´ebrica, obtemos a igualdade: ˆ θT−θ=Idn×n⊗ZT 0 XtX∗ tdt−1 ·Σ1 2⊗Id2n×2n· ZT 0 XtdWn+1,t ZT 0 XtdWn+2,t . . . ZT 0 XtdW2n,t ,(1.9)
26 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO isto ´e ˆ An:2n−An:2n= Σ1 2ZT 0 Xtd(Wn+1,t Wn+2,t ···W2n,t)∗·ZT 0 XtX∗ tdt−1 . A partir de qualquer destas igualdades, ´e imediato verificar que o estimador ˆ θT´e centrado, aplicando as propriedades dos integrais de Itˆo. A consistˆencia e a eficiˆencia de ˆ θTn˜ao fazem intervir a n˜ao invertibilidade da matriz B, pelo que advˆem da teoria geral e podemos considerar que est˜ao demonstradas em [9]. No entanto, a obten¸c˜ao da normalidade assint´otica do EMV n˜ao ´e imediata, devido `a singularidade da matriz de difus˜ao, B. Por esta raz˜ao, o modelo (1.3) – (1.4) fica logo exclu´ıdo na aplica¸c˜ao de resultados e t´ecnicas explorados em [2] (ver [Teorema 4.6.2]) e [83], que se basearam em condi¸c˜oes de regularidade desta matriz. Uma poss´ıvel abordagem `a investiga¸c˜ao da normalidade assint´otica do estimador dado em (1.8), diferente das exploradas em [2] e [83], ´e apresentada na sec¸c˜ao seguinte. Obviamente, no que respeita `a estima¸c˜ao das matrizes KeC, propriamente ditas, uma transforma¸c˜ao linear ´e respons´avel pela obten¸c˜ao dos resultados da estima¸c˜ao destas matrizes a partir dos resultados obtidos para ˆ θT: [ˆ Kˆ C] = −M(vec−1(ˆ θ))∗, onde vec−1representa a opera¸c˜ao de constru¸c˜ao de uma matriz de tamanho apropriado, a partir de um vetor, coluna a coluna. 1.4 Propriedade LAN do EMV A propriedade LAN do estimador ˆ θTdado por (1.8) pode ser obtida usando a Transformada de Laplace de uma martingala que constru´ımos a partir das observa¸c˜oes do processo, seguindo o mesmo tipo de abordagem exposto em [23]. Por outras palavras, reportamo-nos ao uso das condi¸c˜oes de Ibragimov–Khasminskii’s no nosso contexto em particular, que asseguram que a propriedade LAN se verifica (ver [43]). De facto, em [23] demonstra-se que, se a Transformada de Laplace de uma martingala que constru´ımos a partir das observa¸c˜oes do processo (que especificaremos mais adiante no texto), se comporta assintoticamente como a exponencial da matriz de Informa¸c˜ao de Fisher, ent˜ao as condi¸c˜oes de Ibragimov–Khasminskii’s verificamse. Essas condi¸c˜oes, por sua vez, garantem a normalidade assint´otica local do estimador
1.4. PROPRIEDADE LAN DO EMV 27 (Propriedade LAN). O principal resultado desta sec¸c˜ao estabelece essa mesma convergˆencia em lei e ´e o seguinte: Teorema 1.4.1 O EMV definido por (1.8) ´e tal que: √T(ˆ θT−θ)−→ LN0,Σ⊗I−1(A), quando T→+∞,onde I(A)´e solu¸c˜ao da equa¸c˜ao de Lyapunov AI+IA∗+B= 0 .(1.10) De modo a provar este teorema, temos de provar primeiro que a seguinte propriedade se verifica para o processo {Xt}. Daqui para a frente, para simplificar a escrita, usaremos a nota¸c˜ao p= 2n. Proposi¸c˜ao 1.4.1 O processo {Xt}´e tal que: ∀α∈Rp,lim T→+∞Eexp −1 2TZT 0 α∗XtX∗ tαdt= exp −1 2α∗I(A)α, onde a matriz I(A)´e a ´unica solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10) . De modo a provar esta propriedade vamos recorrer a dois lemas que estabelecem resultados auxiliares. Lema 1.4.1 Verifica-se a seguinte convergˆencia: ∀α∈Rp,lim T→+∞Eexp −1 2TZT 0 α∗XtX∗ tαdt= exp −1 2 p X j=1 λ′ j(0) , onde λ′ j(0) designa a derivada com respeito a µde λj(µ)em µ= 0, sendo λjum qualquer valor pr´oprio de Λµtal que ℜe(λj(µ)) >0com Λµ= −A B µαα∗A∗ . Note-se que, para pequenos valores de µ, o espetro de Λµcont´em a parte λj(µ) tal que ℜe(λj(µ)) >0 e a parte λi(µ) tal que ℜe(λi(µ)) <0, aproximadamente sp (Λ0) = {−λj(A)}∪ {λj(A)}, pelo que nos interessaremos pela matriz Λµcom µ=1 Te pela sua convergˆencia quando T→+∞.
28 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO Demonstra¸c˜ao do Lema 1.4.1 Seguindo o mesmo tipo de abordagem usada em [55] definimos a Transformada de Laplace: LT(µ, X) = Eexp −µ 2ZT 0 α∗XtX∗ tαdt. Recorrendo a argumentos idˆenticos aos usados em [23] (ver Sec¸c˜ao 3), podemos escrever que: LT(µ, X) = exp −1 2µtr(A)(det Ψ1(T, µ))−1 2,(1.11) onde [Ψ1Ψ2] ´e a solu¸c˜ao de · Ψ1· Ψ2=hΨ1Ψ2iΛµ hΨ1(0) Ψ2(0)i= [Id 0] , com Λµ= Λ0+µH , Λ0= −A B 0A∗ eH= 0 0 αα∗0 . Isto significa que podemos escrever a seguinte decomposi¸c˜ao: Ψ1(T, µ) = [Id 0] GµD(eT λk(µ))G−1 µ Id 0 sendo Gµassintoticamente triangular superior por blocos e D(·) uma matriz diagonal de dimens˜ao 2p×2p. Descrevendo mais detalhadamente esta decomposi¸c˜ao, temos que a matriz D´e da forma D= D10 0D2 , sendo D1=D(eT λj(µ)) tal que ℜe(λj(µ)) >0 e D2=D(eT λi(µ)) tal que ℜe(λi(µ)) <0. Desta decomposi¸c˜ao conclui-se que Ψ1(T, µ) = C1D1C2+C3D2C4, sendo as matrizes C1,C2,C3eC4, resultantes da pr´opria decomposi¸c˜ao, ou seja, dos c´alculos matriciais, n˜ao sendo relevante, para o que se segue, serem determinadas e interessando apenas as suas propriedades. Usando propriedades alg´ebricas do determinante temos que det Ψ1(T, µ) = det (C1D1C2) det Id + (C1D1C2)−1(C3D2C4).
1.4. PROPRIEDADE LAN DO EMV 29 Por outro lado, tr(A) = Pp i=1 λi(0) = Pp j=1 (−λj(0)) ,para os valores pr´oprios λjtais que ℜe(λj(0)) >0. Logo, usando (1.11), vem que: LT(µ, X) = exp −1 2µ p X j=1 −λj(0) det (C1D1C2) det Id + (C1D1C2)−1(C3D2C4)−1 2 = exp 1 2µ p X j=1 λj(0) exp −1 2µ p X j=1 λj(µ) det Id + (C1D1C2)−1(C3D2C4)−1 2 = exp −1 2µ p X j=1 (λj(µ)−λj(0)) det Id + (C1D1C2)−1(C3D2C4)−1 2 e, usando a expans˜ao de Taylor para os valores pr´oprios λj(µ) de Λµtais que ℜe(λj(µ)) >0: λj(µ) = λj(0) + µλ′ j(0) + o(µ2) quando µ=1 T, obtemos que lim T→+∞LT(X) = exp −1 2 p X j=1 λ′j(0) (1 + o(1)) , ficando conclu´ıda a demonstra¸c˜ao. Lema 1.4.2 Usando a nota¸c˜ao do Lema 1.4.1 verifica-se a seguinte igualdade: p X j=1 λ′j(0) = α∗I(A)α , onde a matriz I(A)´e a ´unica solu¸c˜ao de (1.10). Demonstra¸c˜ao do Lema 1.4.2 Definimos P(λ, µ) = det (Λ0+µH −λId). Aplicando o Teorema da Fun¸c˜ao Impl´ıcita `a equa¸c˜ao caracter´ıstica de Λµ, i.e. P(λ, µ) = 0, podemos facilmente calcular a derivada de λcom respeito a µem 0: λ′ µ(0) = −P′ µ(λ, 0) P′ λ(λ, 0) ,(1.12)
36 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO T200 500 1000 2000 m´edia (ˆ k) / desv. pad. (ˆ k) 35.21/0.405 35.22/0.273 35.22/0.189 35.21/0.142 m´edia (ˆc) / desv. pad. (ˆc) 0.556/0.070 0.565/0.050 0.559/0.032 0.569/0.023 Tabela 1.1: Valor m´edio e desvio padr˜ao obtidos a partir das estimativas dos parˆametros k ec, para diferentes valores de T(Exemplo 1). T200 500 1000 2000 valor −ppara ˆ k0.5336 0.4966 0.6691 0.9900 valor −ppara ˆc0.7447 0.4825 0.8844 0.9447 Tabela 1.2: Resultados do teste de normalidade de Kolmogorov-Smirnov aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de T, para testar a distribui¸c˜ao assint´otica marginal conforme estabelecida no Teorema 1.4.1 (Exemplo 1). O coeficiente de correla¸c˜ao de Pearson entre ˆ ke ˆccom base nos valores obtidos na simula¸c˜ao ´e -0.0723 para T= 2000 s. O teste `a significˆancia do coeficiente de correla¸c˜ao n˜ao apresenta evidˆencias para rejei¸c˜ao da hip´otese de n˜ao correla¸c˜ao entre ˆ ke ˆcpara um n´ıvel de significˆancia de 1%. Os valores de prova para os diferentes valores de Tconstam na Tabela 1.3. T200 500 1000 2000 valor −p0.3101 0.3878 0.4895 0.8627 Tabela 1.3: Resultados do teste de n˜ao correla¸c˜ao para ˆ ke ˆc(Exemplo 1). A Figura 1.1 mostra as regi˜oes de confian¸ca para ˆ ke ˆcobtidas para diferentes valores de T obtidas usando o software R. Podemos observar como as estimativas dos parˆametros kec v˜ao ficando cada vez mais concentradas perto dos verdadeiros valores de kec, conforme T aumenta. N˜ao ´e de estranhar, atendendo a que a variabilidade das estimativas, de acordo com o Teorema 1.4.1, converge quando T→+∞e, quando Taumenta, a variabilidade das estimativas dever´a diminuir. As Figuras 1.2 e 1.3 mostram-nos os histogramas de ˆ ke ˆc obtidos para diferentes valores de Te as curvas a ponteado representam a curva gaussiana obtida a partir do Teorema 1.4.1 enquanto que as curvas a cheio representam o ajustamento
1.6. APLICAC¸ ˜ OES 37 de uma curva gaussiana aos dados. 34.5 35.0 35.5 36.0 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 estimativas de k (kN/m) estimativas de c (kNs/m) (a) 34.5 35.0 35.5 36.0 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 estimativas de k (kN/m) estimativas de c (kNs/m) (b) 34.5 35.0 35.5 36.0 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 estimativas de k (kN/m) estimativas de c (kNs/m) (c) 34.5 35.0 35.5 36.0 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 estimativas de k (kN/m) estimativas de c (kNs/m) (d) Figura 1.1: Regi˜oes de confian¸ca para ˆ ke ˆcobtidas para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros.
38 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO estimativas de k (kN/m) frequência 34.0 34.5 35.0 35.5 36.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 (a) estimativas de k (kN/m) frequência 34.5 35.0 35.5 36.0 0.0 0.5 1.0 1.5 (b) estimativas de k (kN/m) frequência 34.6 34.8 35.0 35.2 35.4 35.6 35.8 0.0 0.5 1.0 1.5 2.0 (c) estimativa de k (kN/m) frequência 34.8 35.0 35.2 35.4 35.6 0.0 0.5 1.0 1.5 2.0 2.5 3.0 (d) Figura 1.2: Histogramas para ˆ kobtidos para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva gaussiana aos dados.
1.6. APLICAC¸ ˜ OES 39 estimativas de c (kNs/m) frequência 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0 5 10 15 20 (a) estimativas de c (kNs/m) frequência 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0 2 4 6 8 10 (b) estimativas de c (kNs/m) frequência 0.50 0.55 0.60 0 5 10 15 (c) estimativas de c (kNs/m) frequência 0.4 0.5 0.6 0.7 0.8 0123456 (d) Figura 1.3: Histogramas para ˆcobtidos para diferentes valores de T(Exemplo 1): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva gaussiana aos dados.
40 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO A curva gaussiana a ponteado exibida nas Figuras 1.2 e 1.3, assim como os valores de prova do teste de Kolmogorov-Smirnov que constam na Tabela 1.2, usam os valores dos parˆametros da distribui¸c˜ao marginal obtida a partir do Teorema 1.4.1. Como podemos ver, as estimativas dos parˆametros kecapresentam um enviesamento consider´avel enquanto o tempo Tn˜ao ´e suficientemente grande, embora o enviesamento quase desapare¸ca no tempo. Para valores de Tmuito elevados podemos dizer que os parˆametros do modelo podem ser estimados com bastante precis˜ao (ver Tabela 1.1) e as distribui¸c˜oes marginais assint´oticas do estimador n˜ao apresentam evidˆencias conducentes `a rejei¸c˜ao da hip´otese de gausseanidade com os parˆametros dados pelo Teorema 1.4.1 (ver Figuras 1.2 e 1.3 e Tabela 1.2). Exemplo 2 Consideramos o modelo (1.3) – (1.4) com k= 4 kN/m,c= 0.5kNs/m, m= 1 ton,σ= 1,∆t= 1 s. Este mesmo exemplo foi considerado em [87], embora com observa¸c˜oes em tempo discreto, pelo que faremos uso dele para permitir um estudo comparativo com os resultados publicados nesse artigo. Os resultados da aplica¸c˜ao do estimador (1.8) est˜ao representados nas figuras 1.4 a 1.6 e nas tabelas 1.4 e 1.5. T200 500 1000 2000 m´edia (ˆ k) / desv. pad. (ˆ k) 3.997/0.148 4.001/0.087 3.999/0.061 3.999/0.044 m´edia (ˆc) / desv. pad. (ˆc) 0.509/0.072 0.497/0.043 0.495/0.032 0.495/0.023 Tabela 1.4: Valor m´edio e desvio padr˜ao calculados a partir de estimativas dos parˆametros kec, obtidos para diferentes valores de T(Exemplo 2). Tal como no exemplo anterior, a Figura 1.4 mostra as regi˜oes de confian¸ca para ˆ ke ˆcobtidas para diferentes valores de T. Podemos observar como as estimativas dos parˆametros kec v˜ao ficando cada vez mais concentradas perto dos verdadeiros valores de kecconforme T aumenta. As Figuras 1.5 e 1.6 mostram-nos os histogramas de ˆ ke ˆcobtidos para diferentes valores de Te as curvas a ponteado representam a curva gaussiana obtida a partir do Teorema 1.4.1 enquanto que as curvas a cheio representam o ajustamento de uma curva gaussiana aos
1.6. APLICAC¸ ˜ OES 41 dados. 3.6 3.8 4.0 4.2 4.4 0.3 0.4 0.5 0.6 0.7 estimativas de k (kN/m) estimativas de c (kNs/m) (a) 3.6 3.8 4.0 4.2 4.4 0.3 0.4 0.5 0.6 0.7 estimativas de k (kN/m) estimativas de c (kNs/m) (b) 3.6 3.8 4.0 4.2 4.4 0.3 0.4 0.5 0.6 0.7 estimativas de k (kN/m) estimativas de c (kNs/m) (c) 3.6 3.8 4.0 4.2 4.4 0.3 0.4 0.5 0.6 0.7 estimativas de k (kN/m) estimativas de c (kNs/m) (d) Figura 1.4: Regi˜oes de confian¸ca para ˆ ke ˆcobtidas para diferentes valores de T(Exemplo 2): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros.
42 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO estimativas de k (kN/m) frequência 3.6 3.8 4.0 4.2 4.4 0.0 0.5 1.0 1.5 2.0 2.5 3.0 (a) estimativas de k (kN/m) frequência 3.8 3.9 4.0 4.1 4.2 0 1 2 3 4 5 (b) estimativas de k (kN/m) frequência 3.8 3.9 4.0 4.1 4.2 0 1 2 3 4 5 6 7 (c) estimativas de k (kN/m) frequência 3.90 3.95 4.00 4.05 4.10 02468 (d) Figura 1.5: Histogramas para ˆ kobtidos para diferentes valores de T(Exemplo 2): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva Gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva Gaussiana aos dados.
1.6. APLICAC¸ ˜ OES 43 estimativas de c (kNs/m) frequência 0.3 0.4 0.5 0.6 0.7 0123456 (a) estimativas de c (kNs/m) frequência 0.35 0.40 0.45 0.50 0.55 0.60 0 2 4 6 8 (b) estimativas de c (kNs/m) frequência 0.40 0.45 0.50 0.55 0.60 0 2 4 6 8 10 12 14 (c) estimativas de c (kNs/m) frequência 0.40 0.45 0.50 0.55 0 5 10 15 20 (d) Figura 1.6: Histogramas para ˆcobtidos para diferentes valores de T(Exemplo 2): (a)T= 200s, (b)T= 500s, (c)T= 1000s, (d)T= 2000s. A linha a ponteado representa a curva Gaussiana obtida a partir do Teorema 1.4.1 enquanto que a linha a cheio representa o ajustamento de uma curva Gaussiana aos dados.
44 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO T200 500 1000 2000 valor −ppara ˆ k0.8500 0.8732 0.9002 0.9398 valor −ppara ˆc0.3724 0.6371 0.7979 0.8730 Tabela 1.5: Resultados do teste de normalidade de Kolmogorov-Smirnov aplicado `as estimativas dos parˆametros kecobtidas para diferentes valores de T, para testar a distribui¸c˜ao assint´otica estabelecida no Teorema 1.4.1 (Exemplo 2). Mais uma vez, o teste `a significˆancia do coeficiente de correla¸c˜ao de Pearson para T= 2000 s n˜ao apresenta evidˆencias para rejei¸c˜ao da hip´otese de n˜ao correla¸c˜ao entre ˆ ke ˆcpara um n´ıvel de significˆancia de 1% (valor −p= 0.7885). ` A semelhan¸ca do exemplo anterior, deparamo-nos com estimativas dos parˆametros kecque apresentam um enviesamento consider´avel enquanto o tempo Tn˜ao ´e suficientemente grande, embora o enviesamento quase desapare¸ca no tempo (ver Figuras 1.5 e 1.6) e, para valores de Tmuito elevados, podemos dizer que os parˆametros do modelo podem ser estimados com bastante precis˜ao (ver tamb´em Tabela 1.4), n˜ao sendo de rejeitar a gaussianidade do estimador com os parˆametros dados pelo Teorema 1.4.1 (ver Tabela 1.5). Salientamos que obtivemos um menor enviesamento nas estimativas dos parˆametros comparativamente a [87], mas no caso em que se consideram as observa¸c˜oes em tempo cont´ınuo. Exemplo 3 Neste ´ultimo exemplo, retomamos os parˆametros do Exemplo 1, ou seja, k= 35.2kN/m,c= 0.57 kNs/m em= 0.933 ton, ∆t= 1 smas vamos alterar o valor de σ, considerando diferentes valores para σ. Nas Tabelas 1.6 e 1.7 est˜ao representadas as m´edias das estimativas dos parˆametros kec, obtidas para diferentes valores de Te de σ. T200 500 1000 2000 σ= 5 35.26 35.25 35.24 35.23 σ= 10 35.22 35.23 35.21 35.22 σ= 20 35.23 35.24 35.23 35.25 Tabela 1.6: Valor m´edio das estimativas do parˆametro kobtido para diferentes valores de T e de σ(Exemplo 3).
1.6. APLICAC¸ ˜ OES 45 T200 500 1000 2000 σ= 5 0.565 0.557 0.556 0.558 σ= 10 0.557 0.556 0.556 0.567 σ= 20 0.566 0.561 0.558 0.559 Tabela 1.7: Valor m´edio das estimativas do parˆametro cobtido para diferentes valores de T e de σ(Exemplo 3). No sentido de ilustrar o comportamento do estimador (1.8) ao longo do tempo, para diferentes valores de σ, nas figuras seguintes (Figuras 1.7 a 1.10) representamos as curvas de m´edias das estimativas dos parˆametros kec, obtido para diferentes valores de σ(Figuras 1.7 e 1.8), bem como os respetivos erros quadr´aticos m´edios (Figuras 1.9 e 1.10). Nessas figuras, as curvas respeitantes a σ= 5 est˜ao representadas a preto, a vermelho as curvas respeitantes a σ= 10 e a azul as curvas respeitantes a σ= 20. 0 500 1000 1500 2000 35.0 35.2 35.4 35.6 35.8 36.0 T (s) média das estimativas de k (kN/m) Figura 1.7: Curvas m´edias das estimativas do parˆametro kobtida para diferentes valores de σ(Exemplo 3).
52 CAP´ ITULO 1. MODELO LINEAR ESTOC ´ ASTICO EDE que traduz o modelo mecˆanico) e a matriz Σ satisfazem a propriedade comutativa da multiplica¸c˜ao ent˜ao uma express˜ao expl´ıcita da matriz de covariˆancia assint´otica do EMV pode ser deduzida facilmente. Para al´em disso, conclu´ımos que se a matriz Σ ´e diagonal a matriz de Informa¸c˜ao de Fisher assint´otica ´e diagonal por blocos. Apresent´amos um estudo baseado em simula¸c˜oes para exemplos de dimens˜ao 2 e 4. Os resultados das simula¸c˜oes ilustram o comportamento assint´otico do EMV e est˜ao de acordo com os resultados te´oricos obtidos.
Cap´ıtulo 2 Modelos lineares com especificidades No Cap´ıtulo 1 abordou-se a quest˜ao da estima¸c˜ao da matriz de deriva, A, do modelo linear estoc´astico de dimens˜ao 2ndefinido por (1.3) – (1.4). Mais concretamente, no contexto da mecˆanica, investigou-se a quest˜ao da estima¸c˜ao das matrizes K(matriz de rigidez do sistema) e C(matriz de amortecimento do sistema), que constituem a matriz de deriva A. Neste cap´ıtulo ´e abordada a quest˜ao da estima¸c˜ao da matriz de deriva, A, do modelo (1.3) – (1.4), mas em modelos que, `a partida, apresentam certas particularidades, como por exemplo, as matrizes KeCserem sim´etricas e/ou diagonais, o que leva a diferen¸cas consider´aveis no problema de estima¸c˜ao de parˆametros. Por este motivo, estes modelos s˜ao tratados num cap´ıtulo distinto. Come¸camos por analisar o modelo (1.3) – (1.4) em dimens˜ao 2ncom as matrizes M,K, Ce Σ1 2diagonais (Sec¸c˜ao 2.1) e seguidamente o caso particular do modelo (1.3) – (1.4) de dimens˜ao 4 em que as matrizes Ke Σ1 2s˜ao diagonais e as matrizes MeCs˜ao sim´etricas (Sec¸c˜ao 2.2). Em ambos os modelos, a t´ecnica utilizada para a obten¸c˜ao do EMV segue de muito perto a que expusemos no Cap´ıtulo 1. No primeiro caso, demonstraremos a propriedade LAN do EMV sem recurso `a Transformada de Laplace visto, neste caso, n˜ao ser vi´avel aplicar esta t´ecnica, pois pretendemos estimar apenas algumas das entradas em cada linha de cada uma das matrizes KeCe n˜ao todas as entradas de cada linha destas matrizes, assim o n´umero de parˆametros a estimar ir´a depender das particularidades do modelo em 53
54 CAP´ ITULO 2. MODELOS LINEARES COM ESPECIFICIDADES estudo. Em alternativa, usaremos as propriedades erg´odicas do processo {Xt}e seguiremos a via cl´assica para dedu¸c˜ao da propriedade LAN. No segundo caso, recorreremos a um estudo baseado em simula¸c˜oes para ilustrar a convergˆencia do estimador para a normal (Sec¸c˜ao 2.4), mas n˜ao deduziremos a propriedade LAN. 2.1 Matrizes de rigidez e de amortecimento diagonais O modelo analisado nesta sec¸c˜ao, como referido, ´e o modelo linear estoc´astico de dimens˜ao 2ndefinido por (1.3) – (1.4) sendo as matrizes M,K,Ce Σ1 2diagonais. Mais precisamente, consideramos, (H1)M=mIdn×ne Σ1 2=σIdn×nem que m > 0 e σ > 0 s˜ao conhecidos e (H2)K=diag (k1, ..., kn) e C=diag (c1, ..., cn) em que ki>0,∀i= 1, ..., n eci>0,∀i= 1, ..., n, (n∈N) s˜ao parˆametros desconhecidos. Tal como no Cap´ıtulo 1, pretende-se estimar os parˆametros k1, ..., knec1, ..., cnno instante de tempo Tcom base nas observa¸c˜oes em tempo cont´ınuo do processo {Xt}, no intervalo de tempo [0, T ]. A solu¸c˜ao do modelo (1.3) – (1.4) ´e o denominado processo de OrnsteinUhlenbeck (processo gaussiano estacion´ario) e, devido ao facto de os parˆametros kieci, i= 1, ..., n serem positivos, sendo m > 0, ´e tamb´em um processo erg´odico (ver [66]). Para a determina¸c˜ao do EMV de hK C i, os passos a seguir s˜ao an´alogos aos do Cap´ıtulo 1. •Constr´oi-se a fun¸c˜ao de log-verosimilhan¸ca, L, que depende dos parˆametros k1, ..., kn ec1, ..., cn. •Seguidamente resolve-se o sistema de equa¸c˜oes normais: ∂L ∂k1= 0 . . . ∂L ∂kn= 0 ∂L ∂c1= 0 . . . ∂L ∂cn= 0 (2.1)
2.1. MATRIZES DE RIGIDEZ E DE AMORTECIMENTO DIAGONAIS 55 para obter o estimador de m´axima verosimilhan¸ca dos parˆametros desconhecidos. ` A semelhan¸ca do que vimos no Cap´ıtulo 1, a fun¸c˜ao de log-verosimilhan¸ca ´e a fun¸c˜ao L(k1, ..., kn, c1, ..., cn) = ln dPX dPW (Xt), onde dPX dPW (Xt) ´e dada por (1.7), sendo M,K,Ce Σ1 2matrizes diagonais. Resolvendo o sistema (2.1) obt´em-se o estimador de m´axima verosimilhan¸ca, que ´e dado por: −ˆ k1 m −ˆc1 m . . . −ˆ kn m −ˆcn m = ZT 0 X2 1,tdt ZT 0 X1,tXn+1,tdt . . . 0 0 ZT 0 X1,tXn+1,tdt ZT 0 X2 n+1,tdt . . . 0 0 ... 0 0 ... ZT 0 X2 n,tdt ZT 0 Xn,tX2n,tdt 0 0 ... ZT 0 Xn,tX2n,tdt ZT 0 X2 2n,tdt −1 · ZT 0 X1,tdXn+1,t ZT 0 Xn+1,tdXn+1,t . . . ZT 0 Xn,tdX2n,t ZT 0 X2n,tdX2n,t .(2.2) O EMV de hK C i,assim como o EMV de A, ´e constru´ıdo usando propriedades da ´algebra matricial, a partir da express˜ao (2.2). O EMV (2.2) ´e centrado (o que ali´as se vˆe muito bem pela express˜ao (2.2), usando as propriedades dos integrais de Itˆo), consistente e eficiente (ver [9]). Pelo facto de o modelo ser erg´odico, ´e sabido que o EMV tem distribui¸c˜ao assint´otica normal (ver [66]). A quest˜ao que se coloca est´a em obter a descri¸c˜ao completa desta distribui¸c˜ao assint´otica, ou seja, os seus momentos de ordem 2, uma vez que o estimador ´e centrado. Usemos a nota¸c˜ao θ2n=−k1 m−c1 m... −kn m−cn m∗. Podemos estabelecer a seguinte propriedade LAN para o EMV (2.2):
56 CAP´ ITULO 2. MODELOS LINEARES COM ESPECIFICIDADES Teorema 2.1.1 Sendo M,K,CeΣ1 2matrizes diagonais, verificando (H1) e (H2), o EMV dado por (2.2) ´e tal que: √T(ˆ θ2nT −θ2n)−→ LN(0, C0), quando T→+∞, onde C0´e uma matriz diagonal de ordem 2n×2ndefinida por C0=diag 2k1c1 m2··· 2kncn m2 2c1 m··· 2cn m. Demonstra¸c˜ao Do modelo (1.3) – (1.4) em dimens˜ao 2n, com as matrizes M,K,Ce Σ1 2 nas condi¸c˜oes do Teorema, recorrendo a alguma ´algebra de matrizes vem que: ZT 0 X1,tdXn+1,t ZT 0 Xn+1,tdXn+1,t . . . ZT 0 Xn,tdX2n,t ZT 0 X2n,tdX2n,t = ZT 0 X2 1,tdt ZT 0 X1,tXn+1,tdt . . . 0 0 ZT 0 X1,tXn+1,tdt ZT 0 X2 n+1,tdt . . . 0 0 ... 0 0 ... ZT 0 X2 n,tdt RT 0Xn,tX2n,tdt 0 0 ... ZT 0 Xn,tX2n,tdt RT 0X2 2n,tdt · −k1 m −c1 m . . . −kn m −cn m +σ ZT 0 X1,tdWn+1,t ZT 0 Xn+1,tdWn+1,t . . . ZT 0 Xn,tdW2n,t ZT 0 X2n,tdW2n,t . Substituindo esta igualdade na express˜ao do estimador (2.2), obt´em-se que: ˆ θ2nT −θ2n=σ ZT 0 X2 1,tdt ZT 0 X1,tXn+1,tdt . . . 0 0 ZT 0 X1,tXn+1,tdt ZT 0 X2 n+1,tdt . . . 0 0 ... 0 0 ... ZT 0 X2 n,tdt ZT 0 Xn,tX2n,tdt 0 0 ... ZT 0 Xn,tX2n,tdt ZT 0 X2 2n,tdt −1
2.1. MATRIZES DE RIGIDEZ E DE AMORTECIMENTO DIAGONAIS 57 · ZT 0 X1,tdWn+1,t ZT 0 Xn+1,tdWn+1,t . . . ZT 0 Xn,tdW2n,t ZT 0 X2n,tdW2n,t .(2.3) Definindo a martingala MT= ZT 0 X1,tdWn+1,t ZT 0 Xn+1,tdWn+1,t . . . ZT 0 Xn,tdW2n,t ZT 0 X2n,tdW2n,t , vem que: ˆ θ2nT −θ2n=σ < MT>−1MT.(2.4) De (2.3) facilmente se conclui a centricidade do EMV e, aplicando as propriedades dos integrais de Itˆo com a decomposi¸c˜ao (2.4), ´e imediato deduzir a rela¸c˜ao entre a matriz de covariˆancia do estimador e a varia¸c˜ao quadr´atica da martingala, MT: Eh√Tˆ θ2nT −θ2n√Tˆ θ2nT −θ2n∗i=Tσ2< MT>−1E(MTM∗ T)< MT>−1 =Tσ2< MT>−1=σ2<1 TMT>−1, visto que, E(MTM∗ T) =< MT>(ver [43], pag. 383). Daqui se deduz que: Eh√Tˆ θ2nT −θ2n√Tˆ θ2nT −θ2n∗i−→ T→+∞σ2I−1(A) (ver [49], pag. 357), sendo a matriz, I(A), a ´unica solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10). Como vimos na Sec¸c˜ao 1.5.2, verificando-se a comutatividade do produto das matrizes M−1K,M−1Ce Σ duas a duas, o que ´e v´alido no caso das matrizes diagonais em estudo, a matriz I(A) ´e a seguinte matriz diagonal por blocos: I(A) = 1 2(M−1K)−1(M−1C)−1Σ 0 01 2(M−1C)−1Σ .
58 CAP´ ITULO 2. MODELOS LINEARES COM ESPECIFICIDADES Nas condi¸c˜oes enunciadas neste Teorema, as propriedades elementares da ´algebra de matrizes diagonais resultam na simplifica¸c˜ao: I(A) = diag m2σ2 2k1c1··· m2σ2 2kncn mσ2 2c1··· mσ2 2cn. Finalmente, as propriedades da inversa de matrizes diagonais permitem concluir a demonstra¸c˜ao. Este resultado constitui, de certa forma, a generaliza¸c˜ao mais natural do problema 2D e, em termos do modelo de mˆecanica estrutural ou de circuitos el´etricos, corresponde a situa¸c˜oes que surgem na pr´atica com alguma frequˆencia. 2.2 Matriz de rigidez diagonal e matriz de amortecimento sim´etrica Nesta sec¸c˜ao, o modelo a analisar consiste no caso particular do modelo (1.3) – (1.4) com matrizes de deriva e de difus˜ao de dimens˜ao 4, tendo, na linguagem da mecˆanica que temos vindo a usar, uma matriz de rigidez diagonal e matrizes de amortecimento e de massa sim´etricas, isto ´e, consideramos o seguinte modelo: dXt= 02×2Id2×2 −M−1K−M−1C Xtdt + 02×202×2 02×2Σ1 2 dWt(2.5) com K= k0 0k , sendo k > 0, C= c1c2 c2c1 , sendo c1, c2>0, M= m1m2 m2m1 , sendo m1, m2>0 e Σ1 2= σ0 0σ , sendo σ > 0.Os parˆametros m1, m2s˜ao supostos conhecidos e os parˆametros k, c1ec2s˜ao os parˆametros a estimar. Com o objetivo de obter o estimador de m´axima verosimilhan¸ca para k,c1ec2,como ´e habitual, come¸camos por explicitar a fun¸c˜ao de log-verosimilhan¸ca, usando (1.7) tal como na sec¸c˜ao anterior.
2.2. MATRIZ DE RIGIDEZ DIAGONAL E MATRIZ DE AMORTECIMENTO SIM´ ETRICA59 Assim, L(k, c1, c2) = σ2k2 m2 1−m2 22m2 1+m2 2X2 1,t −4m1m2X1,tX2,t +m2 1+m2 2X2 2,t +σ2k m2 1−m2 22m2 1c1−2m1m2c2+m2 2c1X2 3,t+ +2 m2c2−2m1m2c1+m2 1c2X3,tX4,t +m2 1c1−2m1m2c2+m2 2c1X2 4,t +σ2k m2 1−m2 22m2 1+m2 2c1−2m1m2c2X2 1,t +2 m2 1+m2 2c2−2m1m2c1X1,tX2,t +m2 1+m2 2c1−2m1m2c2X2 2,t +σ2 m2 1−m2 22h(m1c1−m2c2)2+ (m1c2−m2c1)2X2 3,t +4 (m1c1−m2c2) (m1c2−m2c1)X3,tX4,t +(m1c1−m2c2)2+ (m1c2−m2c1)2X2 4,ti. Resolvendo o sistema: ∂L ∂k = 0 ∂L ∂c1 = 0 ∂L ∂c2 = 0 ,(2.6) obt´em-se o estimador de m´axima verosimilhan¸ca ˆ k ˆc1 ˆc2 =F−1G, (2.7)
60 CAP´ ITULO 2. MODELOS LINEARES COM ESPECIFICIDADES onde a matriz F´e uma matriz sim´etrica de ordem 3 cujos elementos s˜ao dados por: f11 =−1 m2 1−m2 2m2 1+m2 2ZT 0X2 1,t +X2 2,tdt −4m1m2ZT 0 X1,tX2,tdt f12 =−m2 1+m2 2ZT 0 (X1,tX3,t +X2,tX4,t)dt + 2m1m2ZT 0 (X1,tX4,t +X2,tX3,t)dt f13 = 2m1m2ZT 0 (X1,tX3,t +X2,tX4,t)dt −m2 1−m2 2ZT 0 (X1,tX4,t +X2,tX3,t)dt f21 =f12 f22 =−1 m2 1−m2 2m2 1+m2 2ZT 0X2 3,t +X2 4,tdt −4m1m2ZT 0 X3,tX4,tdt f23 =−1 m2 1−m2 2−2m1m2ZT 0X2 3,t +X2 4,tdt + 2 m2 1+m2 2ZT 0 X3,tX4,tdt f31 =−2m2 1+m2 2ZT 0 (X1,tX4,t +X2,tX3,t)dt + 4m1m2ZT 0 (X1,tX3,t +X2,tX4,t)dt f32 =f23 f33 =f22 eG´e uma matriz coluna de dimens˜ao 3 cujos elementos s˜ao dados por: g1=m1ZT 0 X1,tdX3,t +ZT 0 X2,tdX4,t−m2ZT 0 X2,tdX3,t +ZT 0 X1,tdX4,t g2=m1ZT 0 X3,tdX3,t +ZT 0 X4,tdX4,t−m2ZT 0 X4,tdX3,t +ZT 0 X3,tdX4,t g3=m1ZT 0 X4,tdX3,t +ZT 0 X3,tdX4,t−m2ZT 0 X3,tdX3,t +ZT 0 X4,tdX4,t. O EMV de hK C ipode ser escrito a partir de (2.7) `a custa de c´alculos matriciais, sendo consistente (ver [9]) e centrado. A centricidade obt´em-se da seguinte igualdade: ˆ K ˆ C = K C +ZT 0 XtX∗ tdt−1ZT 0 Xtd(W3,t W4,t)Σ1 2 ∗ M∗, aplicando as propriedades dos integrais de Itˆo. Uma vez que o processo Xt´e erg´odico fica garantida a convergˆencia para a distribui¸c˜ao gaussiana e, por isso, restar-nos-ia calcular a matriz de covariˆancia assint´otica do estimador. Neste caso, ao contr´ario do que aconteceu na sec¸c˜ao anterior (Sec¸c˜ao 2.1), n˜ao ´e evidente como se pode usar a matriz, I(A), a ´unica
2.3. OUTROS CASOS DE MATRIZES DIAGONAIS E/OU SIM´ ETRICAS 61 solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10), visto que a matriz de covariˆancia do EMV (2.7) ´e uma matriz quadrada de ordem 3, enquanto que a matriz I(A) ´e uma matriz quadrada de ordem 4. Assim, para determinar a matriz de covariˆancia do EMV dado pela express˜ao (2.7) ´e necess´ario calcular a inversa da matriz quadrada de ordem 3 que aparece na express˜ao (2.7), matriz F, e realizar o produto das matrizes, F−1eG, o que permitir´a explicitar o EMV de ˆ k, ˆc1e ˆc2, obtendo uma express˜ao para este EMV que pode ser usada para a determina¸c˜ao da sua matriz de covariˆancia, seguindo agora, um percurso an´alogo ao do Teorema 2.1.1. S˜ao c´alculos que, em termos computacionais, teremos facilidade em implementar mas que, em termos anal´ıticos se tornam fastidiosos. 2.3 Outros casos de matrizes diagonais e/ou sim´etricas Terminamos a nossa ronda pelos casos especiais referindo que, no caso de estarmos a tratar um problema semelhante ao anterior (modelo (2.5)), mas em que a matriz de rigidez, K, ´e diagonal com os elementos da diagonal eventualmente diferentes, ou seja, K= k10 0k2 , sendo k1, k2>0 parˆametros desconhecidos, depois de determinado o respetivo EMV de ˆ k1, ˆ k2, ˆc1e ˆc2, n˜ao ter´ıamos nenhuma dificuldade em obter um teorema equivalente ao Teorema 2.1.1, visto que o processo ´e erg´odico e ´e clara a liga¸c˜ao entre a matriz de covariˆancia assint´otica do EMV e a matriz de Informa¸c˜ao de Fisher. Tamb´em ´e interessante referir que, no caso de estarmos a tratar um problema semelhante ao (2.5), mas em que a matriz de rigidez, K, ´e sim´etrica com K= k1k2 k2k1 , sendo k1, k2>0 parˆametros desconhecidos, ´e poss´ıvel aplicar a t´ecnica da Transformada de Laplace, de forma semelhante ao que foi feito no Cap´ıtulo 1, Teorema 1.4, visto, neste caso, termos uma linha da matriz Acom todas as entradas a estimar, que se repete imediatamente na linha abaixo. Torna-se, no entanto, obviamente necess´ario refazer a etapa de dedu¸c˜ao da express˜ao do EMV. 2.4 Aplica¸c˜ao Na impossibilidade de apresentar um resultado te´orico sobre a distribui¸c˜ao do EMV (2.7), na situa¸c˜ao espec´ıfica que foi descrita, vamos ilustrar, nesta sec¸c˜ao, o comportamento assint´otico do EMV (2.7) com base em resultados obtidos em simula¸c˜oes.
69 Problemas do tipo ”change-point”s˜ao interessantes pelos desafios probabibil´ısticos que apresentam e pelas suas aplica¸c˜oes, em engenharia estes problemas revestem-se de particular interesse ao investigar situa¸c˜oes de degrada¸c˜ao de uma estrutura. Na Parte II desta tese de Doutoramento, motivados por estas quest˜oes pr´aticas que nos vˆem da engenharia, vamos considerar uma estrutura modelada inicialmente pela equa¸c˜ao diferencial estoc´astica dXt= AXtdt +B1 2dWt, mas em que, ap´os um certo per´ıodo de tempo, ocorre uma mudan¸ca na matriz de deriva, A, que caracteriza o modelo. O objetivo do estudo ser´a estimar os parˆametros de rigidez e de amortecimento da estrutura, que intervˆem na matriz de deriva A, estima¸c˜ao essa que ser´a necess´ario fazer antes e depois da mudan¸ca de regime ocorrer, com base nas observa¸c˜oes do processo ao longo de um intervalo de tempo. Admitindo que as observa¸c˜oes s˜ao feitas em instantes de tempo discretos. Antes de mais vamos abordar este problema considerando que se sabe que de facto ocorreu uma mudan¸ca de regime embora se possa desconhecer em que instante ocorreu. Esse estudo ´e feito no Cap´ıtulo 3. Primeiro estudamos o caso do instante em que ocorreu a mudan¸ca de regime ser conhecido ou observado e depois o caso em que esse instante ´e desconhecido ou n˜ao observado. Neste ´ultimo caso, o instante de mudan¸ca de regime passa a ser tamb´em um dos parˆametros a estimar, sendo este o caso de maior interesse pr´atico em engenharia, pela sua maior complexidade e incerteza envolvida. Em modelos em que h´a suspeita de ocorrˆencia de mudan¸ca de regime ´e de extrema importˆancia a realiza¸c˜ao de um teste de hip´oteses para detetar se efetivamente ocorreu uma mudan¸ca de regime. Para n˜ao alongar demasiado a exposi¸c˜ao, o estudo do teste de hip´oteses adequado a esse problema, ´e apresentado num cap´ıtulo distinto (Cap´ıtulo 4).
70
Cap´ıtulo 3 Modelo estoc´astico com mudan¸ca de regime 3.1 Introdu¸c˜ao Na primeira parte desta tese de Doutoramento, designada Parte I, foi apresentado o estudo do problema da estima¸c˜ao estat´ıstica dos parˆametros da matriz de deriva A, do modelo linear estoc´astico em dimens˜ao 2n, dXt=AXtdt +B1 2dWt, considerando as observa¸c˜oes obtidas em tempo cont´ınuo e tendo AeBdeterminadas propriedades. Como foi referido, em engenharia estrutural, esta equa¸c˜ao diferencial estoc´astica tem sido usada para modelar estruturas com comportamento dinˆamico linear sujeitas a a¸c˜oes externas de car´ater aleat´orio. Com base nas observa¸c˜oes do processo, o Cap´ıtulo 1 incidiu sobre a estima¸c˜ao da matriz de rigidez e da matriz de amortecimento que constituem a matriz de deriva A, ou seja, dos parˆametros desconhecidos no modelo estoc´astico. Acontece que, em engenharia, ´e usual surgirem problemas relacionados com as estruturas ao longo da sua vida, que levam a altera¸c˜oes nos parˆametros da matriz de deriva que modela a estrutura. Veja-se por exemplo [74], em que se estuda as vibra¸c˜oes caracter´ısticas de uma ponte e se mostra como altera¸c˜oes nos parˆametros da matriz de deriva que modela a ponte levam a altera¸c˜oes nas suas vibra¸c˜oes caracter´ısticas. Esse estudo refere, em particular, que um dano estrutural leva a uma redu¸c˜ao no coeficiente de rigidez e, como tal, a uma redu¸c˜ao das frequˆencias naturais da estrutura, algo que ´e bem conhecido dos engenheiros. ´ E de salientar a importˆancia da monitoriza¸c˜ao de 71
72 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME estruturas ao longo da sua vida e desde a sua concep¸c˜ao, pois fornece informa¸c˜ao relacionada com o seu estado de conserva¸c˜ao e a sua seguran¸ca. Na segunda parte desta tese de Doutoramento, motivados por estas quest˜oes pr´aticas que nos vˆem da engenharia, vamos considerar uma estrutura modelada inicialmente pela equa¸c˜ao diferencial estoc´astica dXt=AXtdt +B1 2dWtnas mesmas condi¸c˜oes do Cap´ıtulo 1, mas em que, ap´os um certo per´ıodo de tempo, ocorre uma mudan¸ca na matriz de deriva, A, que caracteriza o modelo. Um modelo adequado a uma situa¸c˜ao deste tipo ser´a um modelo com mudan¸ca de regime, em dimens˜ao 2. O objetivo do estudo ser´a o de estimar os parˆametros de rigidez e de amortecimento da estrutura, que est˜ao presentes na matriz de deriva A, estima¸c˜ao essa que ser´a necess´ario fazer antes e depois da mudan¸ca de regime ocorrer, com base nas observa¸c˜oes do processo ao longo de um intervalo de tempo, admitindo que as observa¸c˜oes s˜ao recolhidas em instantes de tempo discretos. Antes de mais, vamos abordar este problema considerando que se sabe que de facto ocorreu uma mudan¸ca de regime. Esse estudo ´e feito no Cap´ıtulo 3. Com este ponto de partida, primeiro estudamos o caso do instante em que ocorreu a mudan¸ca de regime ser conhecido ou observado (Sec¸c˜ao 3.3) e depois o caso em que esse instante ´e desconhecido ou n˜ao observado (Sec¸c˜ao 3.5). Neste ´ultimo caso, o instante de mudan¸ca de regime passa a ser tamb´em um dos parˆametros a estimar, sendo esta situa¸c˜ao a de maior interesse pr´atico em engenharia, pela sua maior complexidade e incerteza envolvida. Consideramos duas abordagens ao problema de estima¸c˜ao dos parˆametros. Com o objetivo de estimar os parˆametros (de rigidez e de amortecimento) que constituem a matriz de deriva A, antes e depois da mudan¸ca de regime ocorrer, sendo o instante de mudan¸ca de regime conhecido ou observado, usamos as probabilidades de transi¸c˜ao dos estados para construir a fun¸c˜ao de verosimilhan¸ca, esta representa a primeira abordagem ao problema de estima¸c˜ao dos parˆametros. No entanto, este caminho, devido `a complexidade da fun¸c˜ao de verosimilhan¸ca, leva apenas a que se obtenha o estimador de m´axima verosimilhan¸ca dos parˆametros por otimiza¸c˜ao num´erica desta fun¸c˜ao. O esfor¸co computacional revelou-se elevado. Por este motivo, abordou-se o mesmo problema, mas considerando numa primeira fase as observa¸c˜oes do processo em tempo cont´ınuo de modo a poder construir a fun¸c˜ao de verosimilhan¸ca usando a derivada de Radon-Nikodym, e ser poss´ıvel explicitar o estimador de m´axima verosimilhan¸ca dos parˆametros, esta representa a segunda abordagem ao problema
3.2. MODELO COM MUDANC¸A REGIME 73 de estima¸c˜ao dos parˆametros. Seguidamente, uma aproxima¸c˜ao ao estimador de m´axima verosimilhan¸ca em tempo discreto ´e obtida por discretiza¸c˜ao do estimador de m´axima verosimilhan¸ca em tempo cont´ınuo. Esta abordagem, permite explicitar o estimador de m´axima verosimilhan¸ca para o nosso problema original, com a vantagem de ser eficiente em termos computacionais. Abordagens deste tipo j´a foram antes utilizadas na literatura (veja-se [21] e [93]), embora em problemas de estima¸c˜ao diferentes do problema aqui tratado. Devido `as vantagens referidas acima, de abordagem ao problema recorrendo a um problema equivalente em tempo cont´ınuo, vantagens essas que pudemos constatar num estudo com base em simula¸c˜oes e que a pr´opria an´alise dos estimadores permite perceber, opt´amos por estudar o problema em que o instante de mudan¸ca de regime ´e desconhecido ou n˜ao observado pela mesma via j´a exposta (Sec¸c˜ao 3.5.2). Neste caso, mostramos, na tese, como ´e poss´ıvel tamb´em estimar o instante de mudan¸ca de regime, para al´em dos restantes parˆametros. Foram realizados estudos com base em simula¸c˜oes com o objetivo de comparar as duas abordagens ao problema de estima¸c˜ao e de verificar a qualidade das estimativas que se obt´em (Sec¸c˜ao 3.4). Os exemplos simulados correspondem a situa¸c˜oes em que ocorre uma dimui¸c˜ao do parˆametro de rigidez, simulando casos de degrada¸c˜ao de uma estrutura. Um estudo sobre problemas do tipo ”change-point”e suas aplica¸c˜oes ´e apresentado em [11]. Em [26] podemos encontrar um resumo geral do problema de dete¸c˜ao de mudan¸ca de regime e de estima¸c˜ao do ”change-point”e em [28] s˜ao discutidos teoremas limite na an´alise destes problemas. Da an´alise destas obras, podemos constatar, que sempre que se tratam problemas do tipo ”change-point”a fase de estima¸c˜ao deve ser precedida da realiza¸c˜ao de um teste de hip´oteses para dete¸c˜ao de ocorrˆencia de mudan¸ca de regime, revestindo-se assim, de particular importˆancia o estudo da realiza¸c˜ao do referido teste de hip´oteses. Este estudo ´e realizado no Cap´ıtulo 4. 3.2 Modelo com mudan¸ca regime Neste cap´ıtulo consideramos uma estrutura simples modelada por um modelo estoc´astico linear assumindo dois regimes distintos, ambos correspondendo `as equa¸c˜oes do oscilador harm´onico simples com um grau de liberdade sujeito a uma for¸ca externa aleat´oria. Mais concretamente, ambos os regimes correspondem ao modelo estrutural sujeito a a¸c˜oes externas
74 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME aleat´orias: m¨x+c˙x+kx =f(t),(3.1) onde xrepresenta o deslocamento, mrepresenta a massa, krepresenta o coeficiente de rigidez, crepresenta o coeficiente de amortecimento e f(t) representa a for¸ca externa aleat´oria que assumimos poder ser representada por um processo de Wiener. Admitimos que este modelo ´e completamente observado e as observa¸c˜oes s˜ao realizadas em tempo discreto. Seguindo as ideias que apresent´amos na Parte I, preferimos reescrever o modelo (3.1) usando a formula¸c˜ao de Itˆo e a equa¸c˜ao diferencial estoc´astica: d X1,t X2,t = 0 1 −k m−c m X1,t X2,t dt + 0 0 0σ dW1,t dW2,t (3.2) onde {X1,t}representa o deslocamento da estrutura, {X2,t}representa a velocidade do movimento, X0representa a condi¸c˜ao inicial assumida como sendo um vetor constante, e {Wt}´e um processo de Wiener estandardizado de dimens˜ao 2. Tal como j´a t´ınhamos adiantado nas p´aginas anteriores, assumimos, numa primeira fase, que a mudan¸ca de regime da estrutura ocorre num instante conhecido ou observado (Sec¸c˜oes 3.3 e 3.4), e posteriormente, consideramos o caso em que a mudan¸ca de regime ocorre num instante desconhecido ou n˜ao observado (Sec¸c˜ao 3.5). Em ambos os casos, a mudan¸ca de regime traduz-se numa altera¸c˜ao de um dos parˆametros que figuram na matriz de deriva que caracteriza o modelo, conforme se especifica mais abaixo. De facto, se uma estrutura for danificada, os valores de kou de cpodem ficar significativamente alterados e portanto ´e importante investigar o problema da estima¸c˜ao da matriz de deriva, A, que pode n˜ao ser constante no tempo. Torna-se por isso pertinente investigar o problema da estima¸c˜ao dos parˆametros da EDE que define o processo: dXt=AtXtdt +B1 2dWt(3.3) admitindo que a matriz Atpode ser constru´ıda com diferentes parˆametros ao longo do tempo, ou, equivalentemente, estudar o problema de estimar as frequˆencias de vibra¸c˜ao da estrutura antes e depois do instante de mudan¸ca de regime. Na sec¸c˜ao seguinte vamos apresentar o modelo que nos interessa estudar e rever algumas caracter´ısticas da solu¸c˜ao da equa¸c˜ao que o representa. Na Sec¸c˜ao 3.3, o problema de
3.2. MODELO COM MUDANC¸A REGIME 75 estima¸c˜ao dos parˆametros ´e abordado seguindo duas abordagens distintas, que submetemos a uma compara¸c˜ao. Veremos que conseguimos explicitar o EMV dos parˆametros seguindo a segunda abordagem, tanto no caso do instante de mudan¸ca de regime ser conhecido ou observado, como no caso em que o instante de mudan¸ca de regime ´e desconhecido ou n˜ao observado. Neste ´ultimo caso, o instante de mudan¸ca de regime tamb´em ´e estimado. 3.2.1 Modelo linear por tro¸cos Neste trabalho consideramos a estrutura modelada pela EDE definida por (3.3) assumindo dois regimes. Assumimos, numa primeira fase, que a mudan¸ca de regime ocorre num instante conhecido ue que a mudan¸ca de regime afeta a matriz de deriva que caracteriza o modelo. Os dois regimes diferem no valor do coeficiente de rigidez que caracteriza a estrutura, isto ´e, a matriz de deriva que caracteriza o modelo (3.3) ´e definida por: At= 0 1 −kt m−c m ,com kt= k1, t ≤u k2, t > u .(3.4) Isto significa que, ap´os o intervalo de tempo [0, u], o coeficiente de rigidez ktmuda de um valor k1>0 para um valor k2>0. Obviamente assumimos que c, m, σ > 0. O modelo (3.3) – (3.4) inclui o cen´ario de uma estrutura que sofreu uma diminui¸c˜ao de rigidez ap´os um determinado per´ıodo de tempo [0, u] (veja-se [74] e [92], por exemplo). 3.2.2 Solu¸c˜ao da EDE linear por tro¸cos (3.3) – (3.4) Facilmente se obt´em que a solu¸c˜ao da EDE (3.3) – (3.4) ´e dada por: Xt= eA1tX0+Zt 0 e−A1(t−s)B1 2dWs, t ≤u eA2(t−u)Xu+Zt u e−A2(t−s)B1 2dWs, t > u (3.5) onde Ai= 0 1 −ki m−c m ,para i= 1,2. O processo {Xt}´e um processo de difus˜ao n˜ao homog´eneo e, em cada regime, ´e um processo de Markov Gaussiano.
76 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME 3.3 Problema de estima¸c˜ao dos parˆametros Considerando a EDE definida por (3.3) – (3.4), o objetivo ´e estimar a matriz de deriva admitindo que os parˆametros meσs˜ao considerados conhecidos. O instante em que ocorre a mudan¸ca de regime, u, como j´a referimos antes, ser´a primeiramente tamb´em assumido como sendo conhecido ou observado e s´o na Sec¸c˜ao 3.5 retiraremos esta suposi¸c˜ao. Sendo assim, o nosso objetivo ´e estimar os parˆametros k1,k2ec, baseando-nos nas observa¸c˜oes do processo {Xt}em tempo discreto obtidas em instantes igualmente espa¸cados, 0,∆t, 2∆t, ..., T ∆t, realizadas no intervalo de tempo [0, T]. Nesta sec¸c˜ao descrevemos como obter as estimativas de m´axima verosimilhan¸ca dos parˆametros da EDE (3.3) – (3.4). Dito de outra forma, veremos como obter as estimativas de m´axima verosimilhan¸ca dos coeficientes de rigidez e amortecimento de uma estrutura antes e depois da mudan¸ca de regime ocorrer e, como consequˆencia, como estimar mudan¸cas nas frequˆencias naturais do oscilador. Consideramos duas abordagens. A primeira abordagem usa as probabilidades de transi¸c˜ao dos estados para construir a fun¸c˜ao de verosimilhan¸ca (Sec¸c˜ao 3.3.1). Neste caso, dada a complexidade desta fun¸c˜ao, n˜ao ´e poss´ıvel obter uma forma expl´ıcita para o estimador de m´axima verosimilhan¸ca, pelo que este ´e obtido por maximiza¸c˜ao num´erica da fun¸c˜ao de verosimilhan¸ca. A segunda abordagem que usaremos (Sec¸c˜ao 3.3.2) consiste em recorrer `a derivada de Radon-Nikodym para construir a fun¸c˜ao de verosimilhan¸ca para um problema de estima¸c˜ao idˆentico ao formulado atr´as mas em que as observa¸c˜oes do processo s˜ao consideradas em tempo cont´ınuo, passando-se posteriormente `a discretiza¸c˜ao do estimador assim obtido, de modo a poder ser usado com observa¸c˜oes em tempo discreto. Como o instante u´e assumido inicialmente como sendo conhecido, o estimador de m´axima verosimilhan¸ca considerando as observa¸c˜oes em tempo cont´ınuo ´e obtido como no Cap´ıtulo 1 para o caso particular do modelo de dimens˜ao 2 (ver Sec¸c˜ao 1.3.1). O estimador de m´axima verosimilhan¸ca ´e ent˜ao discretizado nos instantes que correspondem `as observa¸c˜oes. Esta segunda abordagem segue, em tra¸cos gerais, o ponto de vista exposto em [20].
3.3. PROBLEMA DE ESTIMAC¸ ˜ AO DOS PAR ˆ AMETROS 77 3.3.1 Abordagem 1: EMV aproximado numericamente Conforme referido acima, esta abordagem usa as probabilidades de transi¸c˜ao dos estados para construir a fun¸c˜ao de verosimilhan¸ca e o estimador de m´axima verosimilhan¸ca ´e obtido numericamente por um procedimento computacional de maximiza¸c˜ao da fun¸c˜ao de logverosimilhan¸ca. Mais concretamente, partindo das observa¸c˜oes em tempo discreto, vamos usar as probabilidades de transi¸c˜ao para construir a fun¸c˜ao de log-verosimilhan¸ca, L, em tempo discreto. Sendo nuo ´ındice da ´ultima observa¸c˜ao anterior ao instante de mudan¸ca de regime, a fun¸c˜ao de log-verosimilhan¸ca ´e dada por: L(k1, k2, c) = nu X k=1 ln p1Xtk|Xtk−1(3.6) + ln p3Xtnu+1 |Xtnu+ n X k=nu+2 ln p2Xtk|Xtk−1, onde p1,p2ep3representam as probabilidades de transi¸c˜ao: •Para i= 1,2, ln piXtk|Xtk−1 =−ln(2π)−1 2ln |Ci| −1 2Xtk−eA(i)(tk−tk−1)Xtk−1∗C−1 iXtk−eA(i)(tk−tk−1)Xtk−1,(3.7) onde Ci= tk R tk−1 eA(i)(tk−r)BeA∗ (i)(tk−r)dr •Para i= 3, ln p3Xtnu+1 |Xtnu =−ln(2π)−1 2ln |C3| −1 2Xtnu+1 −eA(3)(tnu+1−tnu)Xtnu ∗C−1 3Xtnu+1 −eA(3)(tnu+1−tnu)Xtnu ,(3.8) onde eA(3)(tnu+1−tnu)=eA(2)(tnu+1−tnu)eA(1)(tnu−tnu−1)
84 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME 34.5 35.0 35.5 36.0 0.40 0.45 0.50 0.55 0.60 0.65 0.70 estimativas de k1 (kN/m) estimativas de c (kNs/m) (a) 28.5 29.0 29.5 30.0 30.5 31.0 31.5 0.50 0.55 0.60 estimativas de k2 (kN/m) estimativas de c (kNs/m) (b) Figura 3.3: Regi˜oes de confian¸ca para Ehˆ ktieE[ˆc] para um n´ıvel de confian¸ca de 95% para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer (Exemplo 1): (a)T= 200s, (b)T= 600s. O ponto m´edio da nuvem ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. Podemos constatar da an´alise da Figura 3.3 que, depois da mudan¸ca de regime ocorrer, em particular para T= 600 s, as estimativas do parˆametro ktainda se encontram bastantes afastadas do seu verdadeiro valor e, por este motivo, o verdadeiro valor dos parˆametros n˜ao se encontra representado, caindo fora da elipse. 3.4.2 Resultados obtidos na abordagem 2 Os resultados obtidos atrav´es da abordagem 2 s˜ao apresentados na tabela e gr´aficos seguintes (Tabela 3.2 e Figuras 3.4 a 3.6). A Tabela 3.2 mostra as m´edias das estimativas dos parˆametros obtidas antes e depois da mudan¸ca de regime ocorrer. Tal como na abordagem anterior, antes da mudan¸ca de regime ocorrer, as estimativas s˜ao as que correspondem a um modelo linear.
3.4. APLICAC¸ ˜ OES 85 T(s) 200 400 600 ˆ k1(kN/m)35.35 35.24 35.24 ˆ k2(kN/m)- 28.47 28.32 ˆc(kNs/m) 0.55 0.55 0.56 Tabela 3.2: (Exemplo 1) Valor m´edio calculado a partir das estimativas dos parˆametros kte cobtidos para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer. Da an´alise da Tabela 3.2 facilmente se percebe que as estimativas obtidas s˜ao quase centradas nos verdadeiros valores, quer antes quer depois da mudan¸ca de regime ocorrer. A Figura 3.4 cont´em os histogramas das estimativas de ktem instantes antes e depois da mudan¸ca de regime ocorrer. Aos histogramas das estimativas de k1ek2ajusta-se uma distribui¸c˜ao normal (p=0.797 e p=0.687 no teste de Shapiro-Wilk, respetivamente) de notar que as estimativas apresentam coeficientes de varia¸c˜ao bastante pequenos. estimativas de k1 (kN/m) frequência 34.0 34.5 35.0 35.5 36.0 36.5 37.0 0.0 0.2 0.4 0.6 0.8 1.0 (a) estimativas de k2 (kN/m) frequência 27.5 28.0 28.5 29.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 (b) Figura 3.4: (Exemplo 1) Histogramas de ˆ ktpara diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer, e curva gaussiana ajustada aos dados: (a)T= 200s, (b)T= 600s.
86 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME A Figura 3.5 cont´em os histogramas das estimativas de cem instantes antes e depois da mudan¸ca de regime ocorrer. Nota-se algum enviesamento nas estimativas, no entanto, ainda se considera que existe um bom ajuste `a distribui¸c˜ao normal tendo por parˆametros o verdadeiro valor de c (p=0.555) e o coeficiente de varia¸c˜ao das estimativas ´e de 1%. estimativas de c (kNs/m) frequência 0.4 0.5 0.6 0.7 0.8 0 1 2 3 4 5 (a) estimativas de c (kNs/m) frequência 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0 2 4 6 8 10 (b) Figura 3.5: (Exemplo 1) Histogramas de ˆcpara diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer, e curva gaussiana ajustada aos dados: (a)T= 200s, (b)T= 600s. A Figura 3.6 mostra as regi˜oes de confian¸ca para Ehˆ ktieE[ˆc] para um n´ıvel de confian¸ca de 95% em instantes antes e depois da mudan¸ca de regime ocorrer. Como se pode ver, as estimativas est˜ao bastante concentradas na elipse que representa a regi˜ao de confian¸ca e bastante pr´oximas do verdadeiro valor dos parˆametros, representado pelo ponto (k, c) indicado pelo c´ırculo a vermelho.
3.5. INSTANTE DE MUDANC¸A DE REGIME N ˜ AO OBSERVADO 87 34.5 35.0 35.5 36.0 0.4 0.5 0.6 0.7 estimativas de k1 (kN/m) estimativas de c (kNs/m) (a) 27.5 28.0 28.5 29.0 0.45 0.50 0.55 0.60 0.65 estimativas de k2 (kN/m) estimativas de c (kNs/m) (b) Figura 3.6: (Exemplo 1) Regi˜oes de confian¸ca para Ehˆ ktieE[ˆc] para um n´ıvel de confian¸ca de 95% para diferentes valores de T, antes e depois da mudan¸ca de regime ocorrer: (a)T= 200s, (b)T= 600s. O ponto m´edio ´e indicado pelo c´ırculo a azul e o c´ırculo a vermelho indica o verdadeiro valor dos parˆametros. Em suma, ambas as abordagens (1 e 2) tˆem um desempenho bastante razo´avel, no entanto a abordagem 2 mostrou melhores resultados no que respeita `a qualidade das estimativas dos parˆametros. Em termos de custos computacionais, a abordagem 1 usa aproximadamente 98 minutos por simula¸c˜ao, enquanto que a abordagem 2 usa apenas 3 minutos por simula¸c˜ao usando um processador Intel Core i7-3632 QM 2.20 GHz. A abordagem 2 tamb´em mostrou ser menos exigente em termos de complexidade do algoritmo de programa¸c˜ao. Assim sendo ser´a natural dar preferˆencia `a abordagem 2 sobre a abordagem 1. 3.5 Instante de mudan¸ca de regime n˜ao observado A hip´otese segundo a qual o instante em que a mudan¸ca de regime ocorreu, u, ´e conhecido a priori ou observado ´e bastante restritiva na pr´atica. Come¸c´amos por considerar o cen´ario mais simples (Sec¸c˜ao 3.3) e, como veremos mais `a frente, as abordagens anteriores podem
88 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME lidar com o facto de que o instante de tempo upossa ser desconhecido, introduzindo algumas altera¸c˜oes que iremos explorar nesta sec¸c˜ao. 3.5.1 Estima¸c˜ao do instante de mudan¸ca de regime Nesta sec¸c˜ao, assumimos que o instante no qual a mudan¸ca de regime ocorre ´e desconhecido ou n˜ao observado. Por simplicidade de apresenta¸c˜ao, apenas, assumimos que a mudan¸ca de regime ocorre num dos instantes 0,∆t, 2∆t, ..., T ∆t,para os quais o processo ´e observado. N˜ao sendo esse o caso, haver´a pequenas adapta¸c˜oes a fazer, no que se segue. A ideia apresentada em [21] e [93] ´e a de que, aplicando o m´etodo da m´axima verosimilhan¸ca por etapas, podemos estimar simultaneamente o instante em que ocorreu a mudan¸ca de regime, u, bem como os coeficientes desconhecidos da matriz de deriva do modelo, no nosso caso k1,k2ec. Mais precisamente, considerando a fun¸c˜ao de log-verosimilhan¸ca dependendo de k1,k2,ce u, como usual, primeiro calculamos a derivada de L(k1, k2, c, u) em ordem a k1,k2ecapenas e resolvemos o seguinte sistema de equa¸c˜oes: ∂L ∂k1 = 0, ∂L ∂k2 = 0, ∂L ∂c = 0,(3.11) para encontrar ˆ k1(u), ˆ k2(u), ˆc(u), isto ´e, os estimadores de m´axima verosimilhan¸ca para k1, k2ec, sendo estes dependentes de u. Somente depois de conclu´ıdo este passo, encontramos a estimativa da m´axima verosimilhan¸ca do instante de mudan¸ca de regime u, dada por: ˆu=argmaxL ˆ k1(u),ˆ k2(u),ˆc(u).(3.12) As estimativas finais de m´axima verosimilhan¸ca para k1,k2ec, no instante de mudan¸ca de regime, s˜ao finalmente obtidas substituindo em ˆ k1(u), ˆ k2(u), ˆc(u) o parˆametro upela sua estimativa de m´axima verosimilhan¸ca dada por 3.12. 3.5.2 Estimadores de k1,k2,ceu No que se segue recorremos `a segunda abordagem ao problema de estima¸c˜ao, apresentada no Par´agrafo 3.3.2. De acordo com o estudo que expusemos na Sec¸c˜ao 3.5, esta abordagem seria a mais vantajosa.
3.5. INSTANTE DE MUDANC¸A DE REGIME N ˜ AO OBSERVADO 89 Para udesconhecido, considere-se ent˜ao a fun¸c˜ao de log-verosimilhan¸ca sugerida no Par´agrafo 3.3.2: L(k1, k2, c, u) = −1 mσ2k1Zu 0 X1,tdX2,t +k2ZT u X1,tdX2,t +cZT 0 X2,tdX2,t −1 2m2σ2Zu 0 (k1X1,t +cX2,t)2dt +ZT u (k2X1,t +cX2,t)2dt(3.13) O estimador de m´axima verosimilhan¸ca para k1,k2ecj´a foi apresentado em (3.10) e ´e definido, para T > u, por: −ˆ k1(u) m −ˆ k2(u) m −ˆc(u) m = Zu 0 X2 1,tdt 0Zu 0 X1,tX2,tdt 0ZT u X2 1,tdt ZT u X1,tX2,tdt Zu 0 X1,tX2,tdt ZT u X1,tX2,tdt ZT 0 X2 2,tdt −1 Zu 0 X1,tdX2,t ZT u X1,tdX2,t ZT 0 X2,tdX2,t . A discretiza¸c˜ao no tempo a aplicar ´e a que foi usada no Par´agrafo 3.3.2. O estimador de m´axima verosimilhan¸ca de u´e (ver (3.12)): ˆu=argmaxL ˆ k1(u),ˆ k2(u),ˆc(u). ´ E necess´ario recalcular as estimativas de k1,k2ec, isto ´e, obter ˆ k1(ˆu), ˆ k2(ˆu), ˆc(ˆu) como as estimativas finais no instante de mudan¸ca de regime. 3.5.3 Aplica¸c˜oes Nesta sec¸c˜ao pretendemos mostrar os resultados obtidos atrav´es da abordagem 2 quando o instante de mudan¸ca de regime ´e desconhecido ou n˜ao observado. Nas 200 simula¸c˜oes da resposta do processo, o instante de mudan¸ca de regime foi simulado usando a distribui¸c˜ao uniforme no intervalo (250 s, 350 s) e o valor simulado foi mantido idˆentico nas 200 trajet´orias simuladas. Novamente, simul´amos as trajet´orias do processo e estim´amos os parˆametros do modelo e o instante desconhecido de mudan¸ca de regime, seguindo o m´etodo descrito no Par´agrafo 3.5.2. Vamos apresentar nesta tese 4 exemplos (Exemplos 2 a 5) que tˆem como ponto de partida o mesmo exemplo apresentado no in´ıcio da Sec¸c˜ao 3.4 (Exemplo 1), com a diferen¸ca do instante de mudan¸ca de regime ser desconhecido ou n˜ao observado. Os trˆes ´ultimos exemplos deste
90 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME par´agrafo dedicado a aplica¸c˜oes apresentam tamb´em uma menor redu¸c˜ao no coeficiente de rigidez, no que se refere `a mudan¸ca de regime, do que o que aconteceu no primeiro exemplo. Os resultados obtidos nos 4 exemplos s˜ao apresentados nas tabelas e gr´aficos seguintes. Exemplo 2 Consideramos o modelo (3.3) – (3.4) com k1= 35.2 (kN/m), k2= 28.16 (kN/m), c= 0.57 (kNs/m), m= 0.933 (ton), σ= 1, T= 600 (s), udesconhecido. Como referido, consideramos um exemplo com os mesmos valores dos parˆametros que o exemplo apresentado no in´ıcio da Sec¸c˜ao 3.4. A Tabela 3.3 resume os resultados obtidos na estima¸c˜ao, apresentando as m´edias das estimativas dos parˆametros quando o valor obtido na simula¸c˜ao do instante de mudan¸ca de regime foi u= 263 s. T(s) 600 ˆu(s) 263.45 ˆ k1(KN/m)35.17 ˆ k2(KN/m)28.24 ˆc(KN.s/m) 0.565 Tabela 3.3: M´edias das estimativas dos parˆametros quando o instante de mudan¸ca de regime ´e desconhecido (u= 263 s). Por consequˆencia, obtemos as seguintes m´edias para as frequˆencias naturais estimadas: sˆ k1 m= 6.14 (Hz) e sˆ k2 m= 5.5 (Hz). As Figuras 3.7 a 3.13 mostram os resultados obtidos para diferentes valores de Tat´e T= 600 s para apenas uma das 200 trajet´orias simuladas. Constatamos que o instante de mudan¸ca de regime, u, ´e facilmente detetado (ver Figura 3.7) e por isso n˜ao ´e surpreendente que seja estimado com bastante precis˜ao. Contudo, repare-se como a ocorrˆencia da mudan¸ca de regime, que diz respeito ao parˆametro k, afeta a estima¸c˜ao do coeficiente de amortecimento c (Figura 3.8). Na Figura 3.9 est´a representada a fun¸c˜ao de log-verosimilhan¸ca para o Exemplo 2 calculada sobre uma das 200 trajet´orias, onde, mais uma vez, pudemos detetar o instante de mudan¸ca de regime, u, que relembramos, ´e o ponto onde ´e atingido o m´aximo desta fun¸c˜ao. As Figuras 3.10 e 3.11 representam as curvas da m´edia das estimativas dos parˆametros de k ec, respetivamente, calculadas sobre as 200 simula¸c˜oes.
3.5. INSTANTE DE MUDANC¸A DE REGIME N ˜ AO OBSERVADO 91 0 100 200 300 400 500 600 26 28 30 32 34 36 T(s) estimativas de k (kN/m) Figura 3.7: (Exemplo 2) Estimativa do coeficiente de rigidez kcalculadas sobre uma trajet´oria ao longo do tempo T(u= 263 s). 0 100 200 300 400 500 600 0.0 0.2 0.4 0.6 0.8 1.0 T(s) estimativas de c (kNs/m) Figura 3.8: (Exemplo 2) Estimativa do coeficiente de amortecimento ccalculadas sobre uma trajet´oria ao longo do tempo T(u= 263 s).
92 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME 0 100 200 300 400 500 600 −1000 0 1000 2000 3000 4000 T (s) L Figura 3.9: Fun¸c˜ao verosimilhan¸ca para o Exemplo 2 calculada sobre uma das 200 trajet´orias (u= 263 s). 0 100 200 300 400 500 600 26 28 30 32 34 36 T(s) média das estimativas de k (kN/m) Figura 3.10: (Exemplo 2) M´edia das estimativas do coeficiente de rigidez calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 263 s).
3.5. INSTANTE DE MUDANC¸A DE REGIME N ˜ AO OBSERVADO 93 0 100 200 300 400 500 600 0.50 0.55 0.60 0.65 0.70 0.75 0.80 T(s) média das estimativas de c (kNs/m) Figura 3.11: (Exemplo 2) M´edia das estimativas do coeficiente de amortecimento calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 263 s). Como esperado, as estimativas dos parˆametros do modelo afastam-se dos verdadeiros valores desses parˆametros na vizinhan¸ca do instante desconhecido de mudan¸ca de regime, mas apenas por um curto per´ıodo de tempo, recuperando a capacidade de estimar os parˆametros do modelo com pequeno erro de estima¸c˜ao, erro esse que desaparece conforme o tempo aumenta. As Figuras 3.12 e 3.13 mostram o RMSE (ra´ız quadrada do erro quadr´atico m´edio) das estimativas dos parˆametros calculado sobre as 200 simula¸c˜oes. Como seria de esperar, a ra´ız quadrada do erro quadr´atico m´edio (RMSE) diminui rapidamente aproximando-se cada vez mais de zero, mantendo-se muito pr´oxima de zero at´e ao fim do intervalo considerado (Figuras 3.12 e 3.13).
100 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME 0 100 200 300 400 500 600 0.50 0.55 0.60 0.65 0.70 0.75 0.80 T(s) média das estimativas de c (kNs/m) Figura 3.19: (Exemplo 5) M´edia das estimativas do coeficiente de amortecimento calculadas sobre as 200 trajet´orias simuladas do processo ao longo do tempo (u= 329 s). Tal como no exemplo anterior, da an´alise da Tabela 3.6 e das Figuras 3.18 e 3.19, constatamos que as estimativas dos parˆametros k1, k2ecainda est˜ao bastante pr´oximas dos verdadeiros valores dos parˆametros a estimar. Quanto `a estima¸c˜ao do instante de mudan¸ca de regime, u, esta mantem-se com uma grande precis˜ao (desvio padr˜ao de ˆuigual a 0.509). Apesar da diferen¸ca ser subtil, comparando com o Exemplo 3, podemos constatar que quanto mais tarde ocorrer a mudan¸ca de regime mais afastada se encontra a estimativa do parˆametro de amortecimento do seu verdadeiro valor. Essencialmente nota-se nas figuras uma altera¸c˜ao brusca nas estimativas do coeficiente de rigidez que acompanha de muito perto a mudan¸ca de regime e uma convergˆencia relativamente suave das estimativas do coeficiente de amortecimento para o seu verdadeiro valor. O parˆametro u´e estimado com erro muito pequeno.
3.6. CONCLUS ˜ OES 101 3.6 Conclus˜oes Este cap´ıtulo da tese prop˜oe duas abordagens que podem ser seguidas quando se quer estimar os parˆametros de um modelo que se apresenta como linear por tro¸cos. O mesmo ´e dizer que investig´amos duas abordagens para estimar as frequˆencias naturais, de um oscilador harm´onico sujeito a a¸c˜oes externas aleat´orias, para o qual ocorreu uma mudan¸ca no coeficiente de rigidez. O estudo de simula¸c˜ao que realiz´amos sobre exemplos pr´oximos da vida real aponta as vantagens da utiliza¸c˜ao da abordagem baseada na discretiza¸c˜ao do estimador correspondendo ao problema em tempo cont´ınuo, mesmo que as observa¸c˜oes sejam recolhidas em tempo discreto. Esta abordagem permite obter estimadores com menor enviesamento, menores erros de estima¸c˜ao e menores custos computacionais. O m´etodo pode lidar com o caso do instante de mudan¸ca de regime ser desconhecido com pequenas adapta¸c˜oes que sugerimos e explor´amos pela via da simula¸c˜ao num´erica. Neste cap´ıtulo tomou-se sempre como certo que ocorreu efetivamente uma mudan¸ca de regime no coeficiente de rigidez da estrutura (podendo o instante de mudan¸ca de regime ser conhecido ou desconhecido). No entanto, na pr´atica, nem sempre temos tantas certezas sobre se ocorreu ou n˜ao uma mudan¸ca de regime. ´ E importante poder-se realizar um teste de hip´oteses para dete¸c˜ao de ocorrˆencia poss´ıvel de mudan¸ca de regime. Por este motivo, o pr´oximo cap´ıtulo da tese (Cap´ıtulo 4) ir´a ser dedicado ao estudo de um teste de hip´oteses para dete¸c˜ao de ocorrˆencia de mudan¸ca de regime.
102 CAP´ ITULO 3. MODELO ESTOC ´ ASTICO COM MUDANC¸A DE REGIME
Cap´ıtulo 4 Teste de mudan¸ca de regime 4.1 Introdu¸c˜ao Como referido, em engenharia, ´e frequente surgirem problemas relacionados com as estruturas que levam a altera¸c˜oes nos parˆametros da matriz de deriva que modela a estrutura. Algo semelhante acontece nos circuitos el´etricos. Em termos de modela¸c˜ao em tempo cont´ınuo, processos de difus˜ao como o processo de Ornstein-Uhlenbeck (OU) tˆem sido, habitualmente escolhidos para modelar o comportamento destas estruturas em que n˜ao ocorreu deterioriza¸c˜ao ou at´e a estrutura entrar num regime de deterioriza¸c˜ao. Apesar da sua importˆancia nas aplica¸c˜oes, n˜ao s´o em engenharia mas tamb´em noutras ´areas como a ´area financeira ou a biol´ogica na literatura, apenas alguns trabalhos tˆem sido dedicados `a dete¸c˜ao de mudan¸cas nos valores dos parˆametros destes processos de difus˜ao. No cap´ıtulo anterior limit´amo-nos a tratar o problema de estima¸c˜ao de parˆametros, assumindo que ocorreu uma mudan¸ca de regime do processo, devido a uma mudan¸ca na matriz de deriva, A, que caracteriza o modelo (3.3). N˜ao pusemos em causa n˜ao saber se essa mudan¸ca teria de facto ocorrido ou n˜ao. Neste cap´ıtulo pretendemos investigar o importante problema de detetar se efetivamente ocorreu uma mudan¸ca nos valores dos parˆametros da matriz de deriva Aque caracteriza o modelo (3.3), num dado intervalo de tempo, quando n˜ao sabemos se essa mudan¸ca ocorreu de facto ou se os parˆametros da matriz Amantiveram um valor constante em todo o intervalo considerado. Em modelos em que h´a suspeita de ocorrˆencia de mudan¸ca de regime ´e de extrema importˆancia a realiza¸c˜ao de um teste de hip´oteses para 103
104 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME detetar se efetivamente ocorreu uma mudan¸ca de regime. Obviamente a realiza¸c˜ao de um tal teste dever´a preceder a fase de estima¸c˜ao, que se far´a conforme vimos no Cap´ıtulo 3. Para uma revis˜ao e actualiza¸c˜ao dos resultados existentes na literatura relativos a testes de hip´oteses de tipo ”change point”e suas aplica¸c˜oes recomendamos livros como por exemplo [11], [26], [28] e [61], bem como os artigos [40], [45], [46], [52], [53], [64], [76] e [77]. Quando se trata o problema de realiza¸c˜ao de testes de hip´oteses de tipo ”change point”podemos usar uma abordagem param´etrica baseada nas fun¸c˜oes de verosimilhan¸ca (ver [40], [45] e [46]) ou uma abordagem n˜ao param´etrica como os testes de tipo ”CUSUM”, em [64], um teste de tipo ”CUSUM”foi usado para testar a mudan¸ca nos parˆametros de processos de difus˜ao. Os trabalhos [53], [76] e [77] s˜ao uma referˆencia te´orica para o estudo destes problemas de realiza¸c˜ao de testes de hip´oteses de tipo ”change point”considerando processos de difus˜ao erg´odicos com observa¸c˜oes em tempo cont´ınuo. Em [53] e [77] ´e proposto um teste de dete¸c˜ao de mudan¸ca nos valores dos parˆametros de deriva de um processo de difus˜ao observado em tempo cont´ınuo, no entanto, a ideia por detr´as da constru¸c˜ao do teste baseiase sobre uma estat´ıstica de teste (estat´ıstica de teste do tipo ”Cramer-von Mises”) que n˜ao segue a abordagem que faremos nas p´aginas seguintes. Contrariamente aos artigos citados anteriormente dispomos de observa¸c˜oes em tempo discreto e n˜ao em tempo cont´ınuo. Salientamos o trabalho [52] onde s˜ao analisadas as estat´ısticas de teste a usar em testes de hip´oteses de tipo ”change point”para dete¸c˜ao de mudan¸ca nos parˆametros de processos de difus˜ao erg´odicos multidimensionais observados em tempo discreto. No presente trabalho, a abordagem adotada para a constru¸c˜ao do teste de hip´oteses ´e inspirada em [32], onde ´e desenvolvida uma t´ecnica para a constru¸c˜ao de um teste de hip´oteses para a dete¸c˜ao de mudan¸ca de regime em modelos autoregressivos em tempo discreto e, sobretudo, na adapta¸c˜ao desse trabalho que foi feita em [34], onde essa mesma t´ecnica foi adotada em modelos observados em tempo cont´ınuo representados por EDEs revers´ıveis para uma m´edia peri´odica. A leitura dessas obras revela ser interessante analisar a abordagem nelas apresentada e adequar essa abordagem ao nosso problema em concreto, pois, conforme vimos no cap´ıtulo anterior, uma abordagem em tempo cont´ınuo envolvendo uma discretiza¸c˜ao mais tardia do procedimento ´e pass´ıvel de fornecer bons resultados, permitindo uma not´oria redu¸c˜ao do esfor¸co computacional. Sendo assim, vamos seguir a mesma ideia na constru¸c˜ao de um teste de hip´oteses para a ocorrˆencia de mudan¸ca de regime, mantendo o problema
4.2. PROBLEMA DE DETEC¸ ˜ AO DE MUDANC¸A DE REGIME 105 inicialmente como se se baseasse em observa¸c˜oes em tempo cont´ınuo para obter a estat´ıstica de teste e s´o depois entramos na fase de discretiza¸c˜ao. A estat´ıstica de teste ´e obtida atrav´es da raz˜ao das verosimilhan¸cas e prova-se que o supremo da estat´ıstica de teste converge em distribui¸c˜ao para uma ponte browniana bidimensional. Esse ´e o resultado te´orico fundamental do Cap´ıtulo 4 (ver Teorema 4.3.1). Alguns autores, como por exemplo [34], s˜ao de opini˜ao que a obten¸c˜ao da convergˆencia da estat´ıstica de teste para a distribui¸c˜ao de uma ponte browniana bidimensional ´e insuficiente para a sua utiliza¸c˜ao em termos pr´aticos, porque consideram ser necess´ario, nesse caso, simular a referida ponte browniana para determinar os valores cr´ıticos do teste. No entanto, n˜ao partilhamos dessa opini˜ao. Encontramos diversos autores (ver [40], [45], [46] e [51], por exemplo) que se dedicaram a estudar m´etodos num´ericos para o c´alculo da distribui¸c˜ao das pontes brownianas e construiram tabelas para a referida distribui¸c˜ao. Na sec¸c˜ao que se segue (Sec¸c˜ao 4.2) apresentamos a formula¸c˜ao do problema de dete¸c˜ao de mudan¸ca de regime e a abordagem adotada para a sua resolu¸c˜ao. Na Sec¸c˜ao 4.3 realizamos a constru¸c˜ao do teste de hip´oteses seguindo a referida abordagem. Na ´ultima sec¸c˜ao deste cap´ıtulo (Sec¸c˜ao 4.4) apresentamos um exemplo de aplica¸c˜ao do teste de hip´oteses referido acima e ilustramos o procedimento a ser seguido. 4.2 Problema de dete¸c˜ao de mudan¸ca de regime Considere-se a equa¸c˜ao diferencial estoc´astica: dXt=AtXtdt +B1 2dWt X0=x0 ,(4.1) onde At= 0 1 −kt m−c m , com kt= k1, t ≤u k2, t > u , e B1 2= 0 0 0σ . Os parˆametros k1, k2ecs˜ao supostos desconhecidos, meσs˜ao conhecidos, x0´e conhecido e {Wt}´e um processo de Wiener de dimens˜ao 2. O problema de detetar se ocorreu uma mudan¸ca nos valores que constituem as entradas da matriz de deriva At´e investigado neste cap´ıtulo. O processo de difus˜ao {Xt}´e supostamente observado em tempo discreto ao longo do intervalo de tempo [0, T], nos instantes 0,∆t, 2∆t, ..., T ∆t, igualmente espa¸cados. Sendo assim, no modelo (4.1),
106 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME o vetor dos parˆametros desconhecidos ´e o vetor θ= (kt, c) e estamos interessados em testar se ocorreu uma mudan¸ca no chamado coeficiente de rigidez, kt, no intervalo de tempo [0, T]. Sendo X1,...,Xnuma sequˆencia de nobserva¸c˜oes do processo de difus˜ao {Xt}nos instantes 0,∆t, 2∆t, ..., T ∆t, pretendemos testar a hip´otese: H0:u=nversus H1:u=nu+1 < n (4.2) A hip´otese nula afirma que n˜ao ocorreu uma mudan¸ca no coeficiente de rigidez no intervalo de tempo [0, T ], enquanto que a hip´otese alternativa afirma que ocorreu uma mudan¸ca no coeficiente de rigidez no instante de tempo nu+1,pertencente ao intervalo de tempo [0, T ]. Tal como no cap´ıtulo anterior (ver Par´agrafo 3.3.1) nurepresenta o instante observado imediatamente antes da mudan¸ca de regime ocorrer. Como referido na introdu¸c˜ao, a abordagem adotada nesta tese para resolver o problema de testar a hip´otese H0contra a hip´otese H1´e nova mas baseia-se nas ideias exposta em [32], no sentido em que temos um contexto distinto mas consistindo na constru¸c˜ao de um teste de raz˜ao de verosimilhan¸cas. A investiga¸c˜ao que vamos apresentar incide sobre dois aspetos: o primeiro prende-se com o facto de, adotando como estat´ıstica de teste uma raz˜ao de verosimilhan¸cas, ser necess´ario explicitar uma express˜ao para essa estat´ıstica de teste ou para o seu logaritmo o que, como veremos mais `a frente, revela ter algumas dificuldades, sobretudo de implementa¸c˜ao; o segundo aspeto prende-se com a necessidade de obter a distribui¸c˜ao assint´otica dessa estat´ıstica de teste ou desse logaritmo quando T−→ +∞. A raz˜ao de verosimilhan¸ca, R´e uma estat´ıstica definida por: R=supk,c l(k, c) supk1,k2,c,u l(k1, k2, c, u), sendo l(k, c) a fun¸c˜ao de verosimilhan¸ca assumindo que n˜ao ocorre mudan¸ca de regime e l(k1, k2, c, u) a fun¸c˜ao de verosimilhan¸ca que adv´em de (3.6). Logo R= supk,c n X k=1 pXtk|Xtk−1 supk1,k2,c,u nu X k=1 p1Xtk|Xtk−1×p3Xtnu+1 |Xtnu× n X k=nu+2 p2Xtk|Xtk−1 onde p´e a densidade de probabilidade de transi¸c˜ao respeitante `a solu¸c˜ao da EDE (4.1), assumindo que n˜ao ocorre mudan¸ca de regime, e p1,p2ep3s˜ao densidades de probabilidade
4.2. PROBLEMA DE DETEC¸ ˜ AO DE MUDANC¸A DE REGIME 107 de transi¸c˜ao que se consideram quando ocorre mudan¸ca de regime e que est˜ao definidas no Par´agrafo 3.3.1 do Cap´ıtulo 3 (express˜oes (3.7) e (3.8)). O logaritmo de Rpode escrever-se, como habitualmente, devido `as estimativas de MV ˆ k, ˆc, ˆ k1,ˆ k2e ˆu: ln R= n X k=1 ln pˆ k,ˆcXtk|Xtk−1 − nu X k=1 ln p1|ˆ k1,ˆ k2,ˆcXtk|Xtk−1−ln p3|ˆ k1,ˆ k2,ˆcXtnu+1 |Xtnu(4.3) − n X k=nu+2 ln p2|ˆ k1,ˆ k2,ˆcXtk|Xtk−1, uma vez que o supremo das fun¸c˜oes de verosimilhan¸ca ´e atingido na estimativa de m´axima verosimilhan¸ca. Para ser poss´ıvel obter uma express˜ao para o logaritmo da raz˜ao de verosimilhan¸ca ´e necess´ario explicitar as probabilidades de transi¸c˜ao p1,p2ep3. Ora, como vimos no cap´ıtulo anterior (Par´agrafo 3.3.1), na pr´atica, estas express˜oes v˜ao ter de ser calculadas numericamente sobre as observa¸c˜oes, `a custa de um esfor¸co computacional muito elevado (mais de 3 horas de processamento para uma trajet´oria, usando um processador Intel Core i7-3632 QM 2.20 GHz). Em suma, devido ao elevado esfor¸co computacional que envolve o c´alculo da express˜ao (4.3), faz sentido procurar uma outra abordagem que permita reduzir esse volume de c´alculo. Conforme vimos no cap´ıtulo anterior, uma abordagem em tempo cont´ınuo envolvendo uma discretiza¸c˜ao mais tardia do procedimento ´e pass´ıvel de fornecer bons resultados, permitindo uma not´oria redu¸c˜ao do esfor¸co computacional. Sendo assim, vamos seguir a mesma ideia na constru¸c˜ao de um teste de hip´oteses para a ocorrˆencia de mudan¸ca de regime, mantendo o problema inicialmente como se se baseasse em observa¸c˜oes em tempo cont´ınuo para obter a estat´ıstica de teste e s´o depois entramos na fase de discretiza¸c˜ao considerando como tendo dispon´ıveis apenas observa¸c˜oes nos instantes 0,∆t, 2∆t, ..., T ∆t. Vamos seguir a linha de racioc´ınio de [34], inicialmente. ´ E o que fazemos na Sec¸c˜ao 4.3. Mais precisamente, no Par´agrafo 4.3.1 vamos determinar a estat´ıstica de teste; no Par´agrafo 4.3.2 a sua distribui¸c˜ao; no Par´agrafo 4.3.3 vamos descrever o procedimento de teste. Na Sec¸c˜ao 4.4 apresentamos exemplos que ilustram a aplica¸c˜ao do referido procedimento. Observa¸c˜ao 4.2.1 ` A semelhan¸ca dos cap´ıtulos anteriores, nos c´alculos matriciais que se
108 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME seguem, θ,ˆ θeθ0ser˜ao identificados como matrizes coluna, sendo θ= −k1 m −k2 m −c m eˆ θ, θ0com o significado habitual. 4.3 Constru¸c˜ao do teste de hip´oteses em tempo cont´ınuo Tomando como ponto de partida o problema em tempo cont´ınuo, o teste de hip´oteses a realizar ´e formulado da seguinte maneira. Pretendemos testar H0:kt=k1versus H1:kt= k1, t ≤u k2, t > u (4.4) sendo u∈[0, T] um instante de mudan¸ca de regime. Como vimos no Cap´ıtulo 1 (Par´agrafo 1.3.1), sob a hip´otese nula, a fun¸c˜ao de verosimilhan¸ca ´e dada por: l(k, c) = exp −1 mσ2 k T Z 0 X1,tdX2,t +c T Z 0 X2,tdX2,t (4.5) −1 2m2σ2ZT 0 (kX1,t +cX2,t)2dt enquanto que, sob a hip´otese alternativa, a fun¸c˜ao de verosimilhan¸ca ´e dada por (3.13) (ver Cap´ıtulo 3), ou seja: l(k1, k2, c, u) = exp −1 mσ2k1Zu 0 X1,tdX2,t +k2ZT u X1,tdX2,t +cZT 0 X2,tdX2,t −1 2m2σ2Zu 0 (k1X1,t +cX2,t)2dt +ZT u (k2X1,t +cX2,t)2dt A raz˜ao de verosimilhan¸cas ´e definida por R=supk,cl(k, c) supk1,k2,c,ul(k1, k2, c, u). 4.3.1 Determina¸c˜ao da estat´ıstica de teste Seguindo de perto a t´ecnica adotada por [34], no que se segue, consideramos que u=sT, com s∈(0,1) ,e vamos estudar a raz˜ao de log-verosimilhan¸ca em fun¸c˜ao de s: ΛT(s) = −2 ln supk,cl(k, c) supk1,k2,c,sT l(k1, k2, c, sT).(4.6)
4.3. CONSTRUC¸ ˜ AO DO TESTE DE HIP ´ OTESES EM TEMPO CONT´ INUO 109 Vamos mostrar como ´e poss´ıvel obter uma representa¸c˜ao expl´ıcita da estat´ıstica de teste (4.6). Atendendo a que o supremo de cada uma das fun¸c˜oes na raz˜ao de verosimilhan¸cas ´e atingido na estimativa de m´axima verosimilhan¸ca temos: ΛT(s) = −2−1 mσ2ˆ kZT 0 X1,tdX2,t + ˆcZT 0 X2,tdX2,t−1 2m2σ2ZT 0ˆ kX1,t + ˆcX2,t2dt +1 mσ2ˆ k1ZsT 0 X1,tdX2,t +ˆ k2ZT sT X1,tdX2,t + ˆcZT 0 X2,tdX2,t(4.7) +1 2m2σ2ZsT 0ˆ k1X1,t + ˆcX2,t2dt +ZT sT ˆ k2X1,t + ˆcX2,t2dt, onde ˆ k, ˆcs˜ao as estimativas de m´axima verosimilhan¸ca sob a hip´otese H0eˆ k1,ˆ k2e ˆcs˜ao as estimativas de m´axima verosimilhan¸ca sob a hip´otese H1. Podemos deduzir o seguinte resultado te´orico sobre a estat´ıstica (ΛT(s))0≤s≤1,que nos servir´a posteriormente para a determina¸c˜ao da sua distribui¸c˜ao em termos assint´oticos (Par´agrafo 4.3.2). Proposi¸c˜ao 4.3.1 A estat´ıstica de teste (ΛT(s))0≤s≤1dada por (4.6) pode ser representada, sob a hip´otese nula, por: ΛT(s) = −R∗ TQ−1 TRT+R∗ sT Q−1 sT RsT + (RT−RsT )∗(QT−QsT )−1(RT−RsT ), onde Qr=Zr 0 XtX∗ tdt eRr=Zr 0 X1,tdW2,t Zr 0 X2,tdW2,t ∗. Demonstra¸c˜ao Consideremos a seguinte decomposi¸c˜ao para a estat´ıstica de teste (4.7): ΛT(s) = −2 [Λ1T+ Λ2T] (4.8) onde Λ1T=−1 mσ2ˆ kZT 0 X1,tdX2,t + ˆcZT 0 X2,tdX2,t−1 2m2σ2ZT 0ˆ kX1,t + ˆcX2,t2dt e Λ2T=1 mσ2ˆ k1ZsT 0 X1,tdX2,t +ˆ k2ZT sT X1,tdX2,t + ˆcZT 0 X2,tdX2,t +1 2m2σ2ZsT 0ˆ k1X1,t + ˆcX2,t2dt +ZT sT ˆ k2X1,t + ˆcX2,t2dt. O nosso objetivo ser´a agora tratar separadamente estas duas parcelas de (4.8) e simplific´a-las.
116 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME Teorema 4.3.1 Sob a hip´otese nula, a estat´ıstica de teste ΛT(s)dada por (4.6) verifica, para quaisquer valores 0< s1< s2<1, a seguinte convergˆencia: sups∈[s1,s2]ΛT(s)−→ Lsups∈[s1,s2]kW(s)−sW (1)k2 s(1 −s), quando T→+∞,onde W´e um movimento browniano de dimens˜ao 2. Isto significa que a distribui¸c˜ao do supremo da estat´ıstica de teste ΛT(s) dada por (4.6) converge para o quadrado da norma de uma ponte browniana de dimens˜ao 2. Como j´a tivemos ocasi˜ao de referir, alguns autores, como por exemplo [34], consideram que resultados deste tipo s˜ao fr´ageis em termos pr´aticos, pois, ao serem aplicados, levam `a necessidade de simular a ponte browniana, o que representaria um problema de alguma complexidade e n˜ao de todo de divulga¸c˜ao corrente, a ser apoiado por software devidamente testado, etc. N˜ao partilhamos dessa opini˜ao. De facto, recentemente, v´arios autores como por exemplo [45], [46] e [51] dedicaram-se a estudar m´etodos num´ericos eficientes para avaliar Psups1≤s≤s2kW(s)−sW (1)k √t(1−t)≥be publicaram tabelas com os respetivos resultados num´ericos. Nomeadamente, em [45], podemos encontrar a seguinte aproxima¸c˜ao ao c´alculo de Psups1≤s≤s2kW(s)−sW (1)k √t(1−t)≥bpara valores de bsuficientemente grandes: P"sups1≤s≤s2kW(s)−sW (1)k pt(1 −t)≥b# =b3e−1 2b2 √2Γ 3 21 21−3 b2log s2(1 −s1) s1(1 −s2)+ 2b−2+Ob−2. Usaremos este resultado e as respetivas tabelas apresentadas em [45] para aproximar a distribui¸c˜ao da estat´ıstica de teste quando tratarmos da aplica¸c˜ao do teste num exemplo, na Sec¸c˜ao 4.4. Quanto `a demonstra¸c˜ao do Teorema 4.3.1, esta obriga-nos a estabelecer algumas proposi¸c˜oes previamente, envolvendo a convergˆencia de algumas das parcelas que intervˆem na estat´ıstica de teste. Como vimos na sec¸c˜ao anterior, a estat´ıstica de teste (4.6), (ΛT(s))0≤s≤1, pode ser representada, sob a hip´otese nula, pela decomposi¸c˜ao (4.14), que ´e precisamente o que ´e dito na
4.3. CONSTRUC¸ ˜ AO DO TESTE DE HIP ´ OTESES EM TEMPO CONT´ INUO 117 Proposi¸c˜ao 4.3.1, ou seja, ΛT(s) = −ZT 0 XtdW2,t∗ZT 0 XtX∗ tdt−1ZT 0 XtdW2,t +ZsT 0 XtdW2,t∗ZsT 0 XtX∗ tdt−1ZsT 0 XtdW2,t +ZT 0 XtdW2,t −ZsT 0 XtdW2,t∗ZT 0 XtX∗ tdt −ZsT 0 XtX∗ tdt−1 ·ZT 0 XtdW2,t −ZsT 0 XtdW2,t =−R∗ TQ−1 TRT+R∗ sT Q−1 sT RsT + (RT−RsT )∗(QT−QsT )−1(RT−RsT ), onde Qr=Zr 0 XtX∗ tdt eRr=Zr 0 X1,tdW2,t Zr 0 X2,tdW2,t ∗,tal como foi descrito na Proposi¸c˜ao 4.3.1. Vamos utilizar esta decomposi¸c˜ao da estat´ıstica de teste e analisar o comportamento assint´otico de cada uma das parcelas desta representa¸c˜ao em termos da sua distribui¸c˜ao assint´otica. Nesse sentido, antes de realizar a demonstra¸c˜ao do Teorema 4.3.1, vamos estabelecer os seguintes resultados auxiliares que permitem analisar o comportamento assint´otico parcela-a-parcela: Proposi¸c˜ao 4.3.2 e Proposi¸c˜ao 4.3.3. Proposi¸c˜ao 4.3.2 Sob a hip´otese H0,tem-se as seguintes convergˆencias: 1. 1 √TZT 0 XtdW2,t −→ LN(0,I(A)) , quando T→+∞e 2. 1 TZT 0 XtX∗ tdt −→ T→+∞I(A), onde I(A)´e solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10) para o caso n= 1 (ver Cap´ıtulo 1, Sec¸c˜ao 1.4). Note-se que, na Proposi¸c˜ao 4.3.2, no ponto 1, estamos a analisar a convergˆencia da martingala 1 √TRT=1 √TZT 0 X1,tdW2,t ZT 0 X2,tdW2,t ∗e no ponto 2 estamos a analisar a convergˆencia da sua varia¸c˜ao quadr´atica 1 TQT=1 TZT 0 XtX∗ tdt. Relembramos que RTeQT aparecem nas parcelas da estat´ıstica de teste (4.6) e por este motivo estamos a realizar este estudo.
118 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME Demonstra¸c˜ao Esta proposi¸c˜ao decorre do Teorema 1.4.1 (ver Cap´ıtulo 1) no caso n= 1. Sen˜ao vejamos. Atendendo a que, sob a hip´otese H0,da express˜ao (1.9), se pode escrever, no caso em que {Xt}´e bidimensional: √Tˆ θT−θ=σ 1 T T Z 0 XtX∗ tdt −1 1 √T T Z 0 XtdW2,t , o Teorema 1.4.1 estabelece a convergˆencia: σ 1 T T Z 0 XtX∗ tdt −1 1 √T T Z 0 XtdW2,t −→ LN0, σI−1(A), onde I(A) ´e solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10) para o caso n= 1 (ver Cap´ıtulo 1, Sec¸c˜ao 1.4). Dado que 1 TZT 0 XtX∗ tdt ´e a varia¸c˜ao quadr´atica da martingala 1 √T T Z 0 XtdW2,t, tem-se que (ver [82]): 1. 1 √TZT 0 XtdW2,t −→ LN(0,I(A)) , quando T→+∞e 2. 1 TZT 0 XtX∗ tdt −→ T→+∞I(A). Podemos tamb´em deduzir a seguinte vers˜ao funcional da proposi¸c˜ao anterior: Proposi¸c˜ao 4.3.3 Sob a hip´otese H0,tem-se que: 1. 1 √TZsT 0 XtdW2,t −→ LN(0, sI(A)) , quando T→+∞e 2. 1 TZsT 0 XtX∗ tdt −→ T→+∞sI(A),
4.3. CONSTRUC¸ ˜ AO DO TESTE DE HIP ´ OTESES EM TEMPO CONT´ INUO 119 com convergˆencia uniforme quase certa no intervalo (0,1) ,onde I(A)´e solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10) (ver Cap´ıtulo 1, Sec¸c˜ao 1.4). Mais uma vez note-se que, na Proposi¸c˜ao 4.3.3, no ponto 1, estamos a analisar a convergˆencia da martingala 1 √TRsT =1 √TZsT 0 X1,tdW2,t ZsT 0 X2,tdW2,t ∗e no ponto 2 estamos a analisar a convergˆencia da sua varia¸c˜ao quadr´atica 1 TQsT =1 TZsT 0 XtX∗ tdt. Relembramos que RsT eQsT aparecem nas parcelas da estat´ıstica de teste (4.6) e por este motivo estamos a realizar este estudo. Demonstra¸c˜ao 1. Dado que 1 √TRsT =1 √TZsT 0 XtdW2,t=√s1 √sT ZsT 0 XtdW2,t´e uma martingala com respeito `a filtra¸c˜ao gerada pelo processo de Wiener {W2,t}no intervalo [0, T ], tem-se que a matriz de covariˆancia de 1 √TRsT ´e: Cov(1 √TRsT ) = s 1 sT ZsT 0 EX2 1,tdt 1 sT ZsT 0 E(X1,tX2,t)dt 1 sT ZsT 0 E(X1,tX2,t)dt 1 sT ZsT 0 EX2 2,tdt =s1 sT ZsT 0 E(XtX∗ t)dt. Para um qualquer valor de s, esta matriz converge, quando T→+∞, para sI(A) (ver [49] pag. 357). O Teorema Limite Central para martingalas (ver [82]) garante o ponto 1na Proposi¸c˜ao. 2. A demonstra¸c˜ao do ponto 2segue na ´ıntegra a demonstra¸c˜ao da Proposi¸c˜ao 3.3 publicada em [34], uma vez que, embora inserida num contexto diferente, a demonstra¸c˜ao da Proposi¸c˜ao 3.3 em [34] requer apenas a convergˆencia estabelecida no ponto 2da Proposi¸c˜ao 4.3.2 desta tese, sendo todos os passos a seguir idˆenticos aos seguidos em [34]. Demonstra¸c˜ao do Teorema 4.3.1 Estamos agora em condi¸c˜oes de proceder `a demonstra¸c˜ao do Teorema 4.3.1 com a ajuda das Proposi¸c˜oes 4.3.1, 4.3.2 e 4.3.3.
120 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME Recorrendo `a Proposi¸c˜ao 4.3.1 relembramos que: ΛT(s) = −R∗ TQ−1 TRT+R∗ sT Q−1 sT RsT + (RT−RsT )∗(QT−QsT )∗(RT−RsT ) Da Proposi¸c˜ao 4.3.2 obtemos a convergˆencia de 1 √TRTe de 1 TQT: •1 √TRT−→ LN(0,I(A)) ,quando T→+∞e •1 TQT−→ T→+∞I(A), e da Proposi¸c˜ao 4.3.3 obtemos a convergˆencia de 1 √TRsT e de 1 TQsT : •1 √TRsT −→ LN(0, sI(A)) ,quando T→+∞e •1 TQsT −→ T→+∞sI(A), onde I(A) ´e solu¸c˜ao da equa¸c˜ao de Lyapunov (1.10) para o caso n= 1 (ver Cap´ıtulo 1, Sec¸c˜ao 1.4). Uma vez que: 1. WT=I(A)−1 2RT´e um movimento browniano cuja matriz de covariˆancia converge para Id2×2quando T→+∞eR∗ TI(A)−1RT=kWTk2; 2. I(A)−1 2RsT ´e uma martingala cuja matriz de covariˆancia converge para sId2×2quando T→+∞eR∗ sT I(A)−1RsT =kWsk2; podemos escrever que: •R∗ TQ−1 TRT=1 √TR∗ T1 TQT−11 √TRT−→ LN(0, Id),quando T→+∞; •R∗ sT Q−1 sT RsT =1 √TR∗ sT 1 TQsT −11 √TRsT −→ LN(0, sId),quando T→+∞; •(RT−RsT )∗(QT−QsT )∗(RT−RsT )−→ LN(0,(1 −s)Id),quando T→+∞; Aplicando o Teorema de Slutsky (ver [85]) tem-se que ΛT(s) = −R∗ TQ−1 TRT+R∗ sT Q−1 sT RsT + (RT−RsT )∗(QT−QsT )∗(RT−RsT ) −→ L−kW1k2+kWsk s+kW1−Wsk2 1−s=kWs−sW1k2 s(1 −s), continuamente no intervalo [s1, s2]⊆[0, T].
4.3. CONSTRUC¸ ˜ AO DO TESTE DE HIP ´ OTESES EM TEMPO CONT´ INUO 121 Do Teorema do Prolongamento Cont´ınuo resulta que: sups∈[s1,s2]ΛT(s)−→ Lsups∈[s1,s2]kWs−sW1k2 s(1 −s), quando T→+∞, ficando conclu´ıda a demonstra¸c˜ao do Teorema. 4.3.3 Descri¸c˜ao do procedimento de teste Uma vez identificada a estat´ıstica de teste e conhecida a sua distribui¸c˜ao assint´otica, nesta sec¸c˜ao, vamos passar a descrever o procedimento do teste de hip´oteses (4.4). No par´agrafo seguinte (Par´agrafo 4.4) ilustraremos a aplica¸c˜ao do referido procedimento com exemplos. Como ´e habitual, o procedimento de implementa¸c˜ao ou aplica¸c˜ao do teste de hip´oteses ´e composto pelas seguintes etapas: 1. C´alculo do valor da estat´ıstica de teste, para o caso concreto; 2. C´alculo do valor cr´ıtico para o n´ıvel de significˆancia pretendido; 3. Aplica¸c˜ao da regra de decis˜ao do teste. Como vimos no Par´agrafo 4.3.1, a estat´ıstica de teste a usar ´e dada por (4.6), e pode ser expressa do seguinte modo: ΛT=−2−1 mσ2ˆ kZT 0 X1,tdX2,t + ˆcZT 0 X2,tdX2,t−1 2m2σ2ZT 0ˆ kX1,t + ˆcX2,t2dt +1 mσ2ˆ k1Zˆu 0 X1,tdX2,t +ˆ k2ZT ˆu X1,tdX2,t + ˆcZT 0 X2,tdX2,t(4.15) +1 2m2σ2Zˆu 0ˆ k1X1,t + ˆcX2,t2dt +ZT ˆuˆ k2X1,t + ˆcX2,t2dt. Etapa 1 O c´alculo do valor da estat´ıstica de teste num caso concreto ´e realizado, obviamente, com base nos valores observados do processo no intervalo [0, T] e faz intervir as estimativas de m´axima verosimilhan¸ca dos parˆametros, calculadas usando a express˜ao (1.3.1) sob a hip´otese H0e a express˜ao (3.10) sob a hip´otese H1. Por uma quest˜ao de simplicidade, nos exemplos de aplica¸c˜ao que se seguem, o EMV, quer sob a hip´otese H0quer sob a hip´otese H1,foi implementado no software R, sendo os integrais que aparecem nas express˜oes envolvidas calculados usando aproxima¸c˜oes num´ericas como vimos no Par´agrafo 3.3.2.
122 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME Etapa 2 O teste de raz˜ao de verosimilhan¸cas ´e um teste de regi˜ao cr´ıtica do tipo {ΛT≥c}. Para um n´ıvel de significˆancia α, escolheremos o maior valor de ctal que P(ΛT≥c)≤αsob a hip´otese H0. Esta etapa integra o momento em que vamos recorrer ao Teorema 4.3.1. O valor cr´ıtico c´e obtido usando a distribui¸c˜ao da ponte browniana, nomeadamente, usando a Tabela 2 na p´agina 438 em [51]. Transcrevemos abaixo alguns dos valores recolhidos dessa tabela: n´ıvel de significˆancia valor cr´ıtico α= 0.01 1.84273 α= 0.05 1.58379 α= 0.1 1.45399 Tabela 4.1: Valores cr´ıticos para diferentes n´ıveis de significˆancia, α, do teste de hip´oteses (4.4). Etapa 3 Nesta ´ultima etapa aplicamos a regra de decis˜ao do teste. Como ´e habitual, num teste de hip´oteses, vamos comparar o valor obtido para a estat´ıstica de teste ΛTsobre os dados com o respetivo valor cr´ıtico ce verificar se o valor da estat´ıstica de teste se encontra ou n˜ao na respetiva regi˜ao cr´ıtica {ΛT≥c}. Caso perten¸ca `a regi˜ao cr´ıtica rejeitamos a hip´otese H0, caso contr´ario n˜ao rejeitamos a hip´otese H0. Como veremos, no par´agrafo seguinte, a implementa¸c˜ao do teste assim descrito, n˜ao oferece dificuldades e, com base em simula¸c˜oes, podemos avaliar a sua efic´acia na dete¸c˜ao de mudan¸ca de regime. 4.4 Aplica¸c˜oes Vamos novamente considerar os modelos apresentados nos Exemplos 2, 3, 4 e 5 do Cap´ıtulo 3, Par´agrafo 3.5.3 mas, para cada um destes exemplos, vamos considerar que n˜ao sabemos se houve ou n˜ao mudan¸ca de regime no intervalo observado [0,600 s].Vamos portanto realizar
4.4. APLICAC¸ ˜ OES 123 o teste de hip´oteses para dete¸c˜ao de mudan¸ca de regime. Relembramos que esta sequˆencia de exemplos representa uma redu¸c˜ao cada vez mais subtil do coeficiente de rigidez do modelo (parˆametro k). Exemplo 1 Consideramos o modelo (4.1) com k1= 35.2 (kN/m), k2= 28.16 (kN/m), c= 0.57 (kNs/m), m= 0.933 (ton), σ= 1, T= 600 (s). O instante de mudan¸ca de regime simulado foi u= 263 (s) para as 200 trajet´orias. A Tabela 4.2 mostra o valor obtido da estat´ıstica de teste sobre uma trajet´oria escolhida ao acaso, e os valores cr´ıticos para diferentes n´ıveis de significˆancia, necess´arios para a tomada de decis˜ao do teste. valor da estat´ıstica de teste 323.37 valor cr´ıtico (α= 0.01) 1.84273 valor cr´ıtico (α= 0.05) 1.58379 valor cr´ıtico (α= 0.1) 1.45399 Tabela 4.2: (Exemplo 1) Valor da estat´ıstica de teste obtido sobre uma trajet´oria e valores cr´ıticos para diferentes n´ıveis de significˆancia (o valor simulado foi u= 263 s). Como o valor da estat´ıstica de teste se encontra claramente na regi˜ao de rejei¸c˜ao para α= 0.01,ou superior, dada esta trajet´oria, rejeitamos a hip´otese H0para qualquer dos n´ıveis de significˆancia e conclu´ımos portanto sobre a existˆencia de uma mudan¸ca de regime. A Figura 4.1 mostra o histograma dos valores obtidos para a estat´ıstica de teste nas 200 simula¸c˜oes realizadas. Apresentamos igualmente, na Figura 4.2, o histograma das estimativas obtidas do instante de mudan¸ca de regime nas 200 trajet´orias geradas.
124 CAP´ ITULO 4. TESTE DE MUDANC¸ A DE REGIME valores da estatística de teste frequência 100 200 300 400 500 600 0.000 0.001 0.002 0.003 0.004 0.005 0.006 Figura 4.1: (Exemplo 1) Valores da estat´ıstica de teste para as 200 trajet´orias geradas na simula¸c˜ao. estimativas do instante de mudança de regime frequência 262.0 262.5 263.0 263.5 264.0 264.5 265.0 0.0 0.2 0.4 0.6 0.8 1.0 Figura 4.2: (Exemplo 1) Estimativas do instante de mudan¸ca de regime para as 200 simula¸c˜oes (o valor simulado do instante de mudan¸ca de regime foi u= 263 s).
4.4. APLICAC¸ ˜ OES 125 Podemos constatar da an´alise da Figura 4.1 que a distribui¸c˜ao da estat´ıstica de teste ´e unimodal com uma ligeira assimetria. Podemos adiantar que, nestas 200 trajet´orias, a decis˜ao do teste de hip´oteses foi sempre a mesma, a de rejei¸c˜ao da hip´otese H0. Exemplo 2 Consideramos o modelo (4.1) com k1= 35.2 (kN/m), k2= 31.68 (kN/m), c= 0.57 (kNs/m), m= 0.933 (ton), σ= 1, T= 600 (s). O instante de mudan¸ca de regime simulado foi u= 285 (s) para as 200 trajet´orias geradas. Relembramos que este exemplo representa um oscilador que, por algum motivo, sofreu uma diminui¸c˜ao de 10% no coeficiente de rigidez num instante de tempo que n˜ao foi observado, ou seja, sofreu uma menor redu¸c˜ao no coeficiente de rigidez do que ocorreu no Exemplo 1 e, portanto, intuitivamente, seria mais dif´ıcil de detetar se houve ou n˜ao mudan¸ca de regime do oscilador. A Tabela 4.3 mostra o valor obtido da estat´ıstica de teste sobre uma das trajet´orias e os valores cr´ıticos para diferentes n´ıveis de significˆancia. valor da estat´ıstica de teste 64.29 valor cr´ıtico (α= 0.01) 1.84273 valor cr´ıtico (α= 0.05) 1.58379 valor cr´ıtico (α= 0.1) 1.45399 Tabela 4.3: (Exemplo 2) Valor da estat´ıstica de teste obtido sobre uma trajet´oria e valores cr´ıticos para diferentes n´ıveis de significˆancia (o valor simulado do instante de mudan¸ca de regime foi u= 285 s). Como o valor da estat´ıstica de teste se encontra na regi˜ao de rejei¸c˜ao para α= 0.01 e valores superiores rejeitamos a hip´otese H0para n´ıveis de significˆancia baixos como 1%. No entanto ´e not´oria a maior aproxima¸c˜ao do valor da estat´ıstica de teste aos seus valores cr´ıticos. As Figuras 4.3 e 4.4 mostram os histogramas dos valores da estat´ıstica de teste e das estimativas do instante de mudan¸ca de regime obtidos nas 200 trajet´orias geradas.