Modelação e Previsão de Velocidade de Ventos

Modelação e Previsão de Velocidade de Ventos Susana Maria Ferreira Pinto Dissertação de Mestrado apresentada à Faculdade de Ciências da Universidade do Porto Mestrado em Engenharia Matemática 2015 Modelação e Previsão de Velocidade de Ventos Susana Maria Ferreira Pinto MSc FCUP ANO 2.º CICLO Modelação e Previsão de Velocidade de Ventos Susana Maria Ferreira Pinto Mestrado em Engenharia Matemática Departamento de Matemática 2015 Orientadora Professora Doutora Margarida Brito, FCUP Todas as correções determinadas pelo júri, e só essas, foram efetuadas. O Presidente do Júri, Porto, ______/______/_________ 4 FCUP Modela¸c˜ao e Previs˜ao de Velocidade de Ventos FCUP 5 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Agradecimentos “Some mathematician, I believe, has said that true pleasure lies not in the discovery of truth, but in the search for it.” — Tolstoi Um especial obrigada ... •Para os meus pais, com quem tenho aprendido a ser quem sou e a quem devo ter chegado at´e aqui. Por toda a paciˆencia e confian¸ca. Por todos os valores e ensinamentos. Por serem as minhas principais fontes de inspira¸c˜ao. •Para a minha Orientadora, Professora Doutora Margarida Brito, por ser minha Professora e conselheira e me transformar numa melhor Matem´atica. Deixo tamb´em um agradecimento a todos os meus professores de Licenciatura em Matem´atica e Mestrado em Engenharia Matem´atica da Faculdade de Ciˆencias da Universidade do Porto, pelo seu contributo no meu percurso acad´emico. •Para aqueles que, de uma forma muito especial, marcaram a minha vida pessoal e universit´aria: Filipe Oliveira, Carla Gon¸calves, Lu´ıs Ba´ıa, Ana Silva, Jos´e Pedro Silva, Filipa Carvalho, J´ulio Silva, Ana Lu´ısa Lopes, Ricardo Cruz, Marisa Reis e Renato Fernandes. Porque todos n´os temos em comum dois aspetos: o gosto pela Matem´atica e uma pan´oplia de aventuras para contar. •Para as amigas de uma vida: Mafalda de Castro, Ana Barros e Joana Freitas. Porque sempre estiveram, est˜ao e estar˜ao l´a para mim. •Para o An´ıbal Couto, que me deu as primeiras luzes sobre o mundo dos seguros. •Para os meus colegas de empresa, por acreditarem em mim e no meu trabalho. 6 FCUP Modela¸c˜ao e Previs˜ao de Velocidade de Ventos FCUP 7 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Resumo Segundo a Organiza¸c˜ao Mundial de Meteorologia, “´e necess´ario lembrar que mesmo as pequenas obstru¸c˜oes podem causar graves altera¸c˜oes na velocidade do vento e desvios na dire¸c˜ao do mesmo”. A determina¸c˜ao precisa da distribui¸c˜ao de probabilidade dos valores de velocidade do vento ´e muito importante na estimativa de potencial energ´etico da velocidade do vento sobre uma regi˜ao. A modela¸c˜ao adequada de ventos e estudo da importˆancia de vari´aveis auxiliares sobre essa modela¸c˜ao (´epoca do ano, existˆencia de precipita¸c˜ao, etc.) pode ser fulcral para entender o que leva, por exemplo, `a baixa rentabilidade de uma vento´ınha e´olica ou at´e ao poss´ıvel descarrilamento de um comboio por combina¸c˜ao de ventos fortes e dire¸c˜ao adversa do mesmo. A previs˜ao com um grau de certeza consider´avel permite a sua utiliza¸c˜ao na meteorologia, aplica¸c˜ao de pr´aticas agr´ıcolas, preven¸c˜ao de ru´ına de companhias seguradoras. Tendo em conta as diferentes aplica¸c˜oes e utilidade deste tipo de estudos, na presente tese s˜ao considerados os processos de modela¸c˜ao e previs˜ao de velocidades m´aximas de vento de uma esta¸c˜ao meteorol´ogica italiana. Tem por base o modelo semi-markoviano para a modela¸c˜ao e faz a sua compara¸c˜ao com o modelo markoviano habitualmente usado. Estuda tamb´em a previs˜ao da referida vari´avel por ´arvores de regress˜ao, redes neuronais, regress˜ao com m´aquinas de suporte vetorial e um modelo linear de s´eries temporais. Finalmente refere um modelo te´orico de modela¸c˜ao de quantidades de ru´ına de uma seguradora — como referˆencia a uma poss´ıvel aplica¸c˜ao do processo modelado. A ideia geral ´e perceber at´e que ponto ´e que um processo real — da´ı o uso de uma base de dados de acesso livre — pode ser explorado com diferentes mecanismos de estudo matem´aticos, abordando ´areas como Data Mining, Estat´ıstica,Teoria de Risco,Modela¸c˜ao Matem´atica,Processos Estoc´asticos eSimula¸c˜ao Computacional. 8 FCUP Modela¸c˜ao e Previs˜ao de Velocidade de Ventos FCUP 9 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Abstract According to World Meteorological Organization, “It must always be remembered that even small obstructions cause serious changes in wind speed and deviations in wind direction”. The precise computation of the probability distribution of wind speed values is very important in the estimation of energy potential of wind speed over a region. The proper modelling of wind and study of the importance of auxiliary variables on this modelling (time of year, the presence of precipitation, and so on...) may be crucial to understand what takes, for instance, low profitability of a wind fan or to a possible derailment of a train by a combination of strong winds and its direction. Forecast with a considerable degree of certainty allows its use in meteorology, application of agricultural practices, prevention of ruin of an insurance company. Taking into account the different applications and utility of such studies, in this thesis are considered processes of modelling and prediction of maximum wind speed of an Italian weather station. It uses the semi-Markovian model for modelling and makes the comparison with the Markovian model commonly used. It also studies the forecast of the variable referred with regression trees, artificial neural networks, regression with support vector machines and a linear time series model. Finally reviews a theoretical model of ruin amounts for an insurance company — reference to a possible application of the modelled process. The general idea is to realize to what extent is that a real process — hence the use of an open access database — can be exploited with different mathematical mechanisms of study, addressing areas such as Data Mining,Statistics,Risk Theory,Mathematical Modelling, Stochastic Processes and Computational Simulation. 16 LISTA DE FIGURAS FCUP 17 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Cap´ıtulo 1 Introdu¸c˜ao A for¸ca da Natureza ´e impar´avel. Lembrando que os eventos clim´aticos a ela associados podem ter consequˆencias nefastas (ou, pelo contr´ario, constitu´ırem um meio para benef´ıcio de atividades humanas), torna-se necess´ario estud´a-los, percebendo como ´e que ´e poss´ıvel model´a-los eprevˆe-los, uma vez que, quando ocorrem, podem ter grandes impactos materiais e pessoais. Para tal, foi feito o estudo de uma base de dados que apresenta diversa informa¸c˜ao hor´aria registada ao longo de dois anos na esta¸c˜ao Meteorol´ogica de Settala (It´alia) e que pode ser encontrada no site http://www.lsi-lastem.com 1. O principal foco ser´a o estudo de velocidades m´aximas de vento. H´a diversas aplica¸c˜oes associadas a este tipo de estudo, como por exemplo, produ¸c˜ao de energia e´olica, defini¸c˜ao de pr´aticas agr´ıcolas, previs˜ao meteorol´ogica, entre outras [Aigner e Gjengedal (2011), Chi et al. (2007)]. Os processos de modela¸c˜ao e previs˜ao de eventos clim´aticos tˆem contribu´ıdo para uma melhoria significativa de estrat´egias de planeamento de atividades frequentemente afetadas pela variabilidade clim´atica, como energia, agricultura e sa´ude [Green et al. (2009)]. A utiliza¸c˜ao dos resultados provenientes da modela¸c˜ao e previs˜ao de eventos clim´aticos obriga ao desenvolvimento de t´ecnicas ou m´etodos que melhorem e aprimorem este tipo de trabalho. Al´em disso, ´e de conhecimento geral que a realidade muitas vezes n˜ao se assemelha `as simplifica¸c˜oes da teoria e que, por isso, as simplifica¸c˜oes te´oricas podem ser demasiadas para conseguir modelar adequadamente um determinado evento. Do mesmo modo, ´e comum ver exemplos de dados muito bem comportados para os quais os modelos funcionam muito bem. E com dados reais, ´e tamb´em esse o comportamento que se deve esperar? Portanto existe, nesta tese, n˜ao s´o uma introdu¸c˜ao de t´ecnicas usadas para os diferentes objetivos tra¸cados anteriormente, mas tamb´em uma aplica¸c˜ao dos mesmos a dados reais. Foi realizada uma separa¸c˜ao f´ısica da informa¸c˜ao, neste documento, em trˆes cap´ıtulos fundamentais: os dois primeiros fazem a contextualiza¸c˜ao te´orica inicialmente pensada para o caso pr´atico, definido no cap´ıtulo seguinte. 1´ultima consulta realizada em 28 de Junho de 2015. 18 CAP´ ITULO 1. INTRODUC¸ ˜ AO FCUP 19 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Cap´ıtulo 2 Modela¸c˜ao “. . . all models are approximations. Essentially, all models are wrong, but some are useful. However, the approximate nature of the model must always be born in mind. . . ” — George E.P. Box As cadeias de Markov [Ross (2014)] s˜ao o processo mais utilizado para modela¸c˜ao dos ventos. Uma das grandes dificuldades que a modela¸c˜ao de ventos evidencia, quando feita desta forma, ´e a falta de flexibilidade que as cadeias de Markov apresentam em rela¸c˜ao aos tempos de espera, o que faz com que as simula¸c˜oes dos ventos sejam pouco fi´eis aos dados originais. Os processos semi-markovianos tˆem sido amplamente utilizados para modelar fen´omenos naturais [Asaduzzaman e MahbubLatif (2013), Barbu et al. (2004), Sansom et al. (2001)] e formam uma classe de processos estoc´asticos que generalizam, em simultˆaneo, cadeias de Markov e processos de renovamento. A principal vantagem da sua utiliza¸c˜ao, quando comparados com os processos Markovianos, ´e o uso de qualquer distribui¸c˜ao para modelar os tempos de espera. 2.1 Breve revis˜ao de processos de Markov Para contextualizar o problema, vai ser feita uma breve revis˜ao de processos estoc´asticos. Um processo estoc´astico ´e uma fam´ılia ou conjunto de vari´aveis aleat´orias definidas num espa¸co de probabilidade - e indexadas no tempo ou no espa¸co. Assuma-se que {X(t)}´e um processo aleat´orio indexado no tempo t∈TeX(t) ´e o estado do processo. Se o conjunto T´e finito ou infinito numer´avel, ent˜ao {X(t)}´e um processo em tempo discreto. Caso contr´ario, ´e um processo em tempo cont´ınuo. Do mesmo modo, se {X(t)}´e definido sob um conjunto cont´avel de estados, ent˜ao o processo ´e definido por uma cadeia, uma vez que ´e discreto nos estados. Caso contr´ario, ´e cont´ınuo nos estados e designado como uma sequˆencia. 20 CAP´ ITULO 2. MODELAC¸ ˜ AO Considerando x(t) como as realiza¸c˜oes de X(t) e assumindo que x(t)∈R, podem definir-se as seguintes fun¸c˜oes: •m´edia como fun¸c˜ao do tempo: µX(t) = mX(t) = E[X(t)] = R∞ −∞ xfX(t)(x(t))dx onde fX(t)representa a fun¸c˜ao densidade de probabilidade de X(t) •autocorrela¸c˜ao de segunda ordem: RX(t1, t2) = E[X(t1)X(t2)] = Z∞ −∞ x1x2fX(t1),X(t2)(x1(t1), x2(t2))dx1dx2 •autocovariˆancia: CX(t1, t2) = Cov(X(t1), X(t2)) = RX(t1, t2)−µX(t1)µX(t2) (observe-se que V ar(X(t)) = CX(t, t)) Sob um ponto de vista de classifica¸c˜ao, e definindo a fun¸c˜ao de distribui¸c˜ao conjunta como FX1(t1),...,Xn(tn)(x1, ..., xn) = P(X(t1)≤x1, ..., X(tn)≤xn), um processo estoc´astico pode ser classificado como: •estacion´ario de ordem nse os momentos de ordem ≤ns˜ao independentes do tempo; •estacion´ario em sentido estrito se, ∀n, ∀{t1, t2, ..., tn},∀τ∈R FX1(t1),...,Xn(tn)(x1, ..., xn) = FX1(t1+τ),...,Xn(tn+τ)(x1, ..., xn); •estacion´ario em sentido lato se µX(t) = E(X(t)) = µXeRX(t+τ, τ) = RX(τ)τ=t1−t2 ou seja, se a m´edia for constante ao longo do tempo e a fun¸c˜ao de autocorrela¸c˜ao s´o depender da diferen¸ca entre t2et1. Note-se que um processo estacion´ario em sentido lato n˜ao implica que o mesmo o seja em sentido estrito, uma vez que a primeira propriedade ´e menos exigente que a segunda. Seja {N(t)}oprocesso de contagem do n´umero de eventos no intervalo (0, t] e N(0) = 0 ≤N(t1)≤N(t2)≤... ≤N(tk)≤... ∀0≤t1≤... ≤tk≤... Lembrando que N(t) = Pn≥tI{Tn≤t},onde I{Tn≤t}=1 se Tn≤t 0 se Tn> t , note-se que {Tn, n ∈N} ≡ {N(t), t ≥0} porque N(t) = n⇔Tn≤t<Tn+1 eN(t)≥n⇔Tn≤t. Esta dualidade ´e muito ´util, porque permite escrever um problema `a custa do outro, caracter´ıstica frequentemente utilizada no estudo deste tipo de processos. 2.2. INTRODUC¸ ˜ AO AOS PROCESSOS SEMI-MARKOVIANOS 21 Diz-se que o processo de contagem, {N(t), t ≥0}tem: •incrementos independentes se N(0), N(t1)−N(0), N(t2)−N(t1), ..., N(tn)−N(tn−1) s˜ao vari´aveis aleat´orias independentes, o que significa que o n´umero de eventos ocorridos at´e ao instante t´e independente do n´umero de eventos que ocorrem entre tet+s, para algum s. •incrementos estacion´arios, se P(N(t+s)−N(t) = k) = P(N(s) = k) para qualquer t≥0. Defini¸c˜ao: Diz-se que um processo {N(t), t ≥0}´e de Poisson se tem incrementos estacion´arios e independentes, N(0) = 0, P(N(h)=1)=λh +o(h) e P(N(h)≥2) = o(h). Proposi¸c˜ao: Se {N(t), t ≥0}´e um processo de Poisson, ent˜ao N(t)∼Poisson(λt), onde λ´e o parˆametro correspondente `a intensidade do processo. Assumindo que {N(t)}segue um processo de Poisson, a distribui¸c˜ao dos tempos entre eventos de uma cadeia de Markov tem distribui¸c˜ao exponencial (a explica¸c˜ao desta afirma¸c˜ao pode ser encontrada no Anexo A). Um processo de Markov ´e um processo estoc´astico caracterizado pela propriedade markoviana, segundo a qual, a probabilidade de qualquer comportamento futuro do processo, quando o seu estado atual ´e conhecido, n˜ao ´e alterada pela existˆencia de conhecimento adicional sobre o seu comportamento passado. Esta informa¸c˜ao ´e traduzida pela seguinte express˜ao: P(X(tk+1)≤xk+1 |X(tk) = xk, ..., X(t0) = x0) = P(X(tt+1)≤xk+1 |X(tk) = xk). 2.2 Introdu¸c˜ao aos processos Semi-markovianos Um processo estoc´astico de Semi-Markov ´e uma generaliza¸c˜ao de um processo de Markov, uma vez que a informa¸c˜ao sobre o tempo de permanˆencia no estado atual passa a ser relevante. Contudo, mant´em-se irrelevante para o comportamento futuro qualquer informa¸c˜ao sobre os estados visitados no passado. Consequentemente, os tempos entre acontecimentos sucessivos deixam de estar restritos `a distribui¸c˜ao exponencial, podendo seguir qualquer distribui¸c˜ao de probabilidade. Seguidamente apresentar-se-˜ao os conceitos fundamentais para a defini¸c˜ao e compreens˜ao de um processo deste tipo, mas uma introdu¸c˜ao mais detalhada pode ser encontrada, por exemplo, em Iosifescu et al. (2013) e Janssen e Manca (2006). Considere-se I={1, ..., M}como um espa¸co finito de estados e (Ω,F,P) como um espa¸co de probabilidade. Sejam Jn: Ω →IeTn: Ω →N 22 CAP´ ITULO 2. MODELAC¸ ˜ AO duas vari´aveis aleat´orias, sendo que Jnrepresenta o estado ocupado pelo sistema imediatamente ap´os ter efetuado a n-´esima transi¸c˜ao e Tnrepresenta o instante de tempo em que ocorreu a n-´esima transi¸c˜ao (n∈N). Xn=Tn−Tn−1´e o tempo que o sistema levou a transitar de Jn−1para Jn. Deste modo, (Xn)n∈Ndesigna uma sequˆencia de tempos de espera e (Tn)n∈Numa sequˆencia de tempos de chegada. Note-se que N(t) = sup{n:Tn≤t} ∀t∈Nrepresenta o n´umero de transi¸c˜oes que ocorreram em (0, t]. Por exemplo, considerando M= 2, uma poss´ıvel realiza¸c˜ao do processo estoc´astico bidimensional (Jn, Xn) seria aquela em que o sistema come¸ca no estado J0= 2, transita para 1 ap´os 2 unidades de tempo (J1= 1, X1= 2) e transita novamente para 2 ap´os 3 unidades de tempo (J2= 2, X2= 3). O par (Jn, Tn) designa o processo Markoviano n˜ao-homog´eneo de renovamento. Note-se que os processos de renovamento fornecem modelos te´oricos de investiga¸c˜ao para a ocorrˆencia de padr˜oes em experiˆencias independentes e repetidas [Pyke (1961)]. De um modo informal, pode escrever-se que o termo renovamento vem do pressuposto b´asico de que, quando o padr˜ao de interesse ocorre pela primeira vez, o processo se repete, no sentido em que a situa¸c˜ao inicial ´e restabelecida. Introduzidos os conceitos b´asicos necess´arios para a compreens˜ao de um modelo semimarkoviano, introduz-se o n´ucleo n˜ao-homog´eneo semi-markoviano associado a este processo, Q= [Qij(s, t)], que representa a probabilidade de chegada a um determinado estado at´e um determinado instante, sabendo qual o estado e instante da transi¸c˜ao anterior, e ´e definido como Qij(s, t) = P(Jn+1 =j, Xn+1 ≤t−s|Jn=i, Tn=s) =P(Jn+1 =j, Tn+1 ≤t|Jn=i, Tn=s) portanto P= [pij(s)] define a matriz de transi¸c˜oes, onde pij(s) = P(Jn+1 =j|Jn=i, Tn=s) devolve a probabilidade do sistema estar no estado jsabendo que no instante sestava no estado i. Note-se ainda que pij(s) = lim t→∞Qij(s, t)i, j ∈I;s, t ∈N, s ≤t. Para al´em do n´ucleo Q, existem outras probabilidades condicionadas relevantes para a defini¸c˜ao de um processo semi-markoviano. Uma fun¸c˜ao particularmente relacionada com a constru¸c˜ao do n´ucleo ´e a fun¸c˜ao de distribui¸c˜ao condicional do tempo de espera em cada estado i, dado o estado subsequente ocupado: Fij(s, t) = P(Xn+1 ≤t−s|Jn=i, Jn+1 =j, Tn=s) 2.2. INTRODUC¸ ˜ AO AOS PROCESSOS SEMI-MARKOVIANOS 23 =P(Tn+1 ≤t|Jn=i, Jn+1 =j, Tn=s) = (Qij (s,t) pij (s)se pij(s)6= 0 1 se pij(s)=0 Note-se que esta fun¸c˜ao pressup˜oe que uma certa transi¸c˜ao vai, de facto, ocorrer, devolvendo apenas a probabilidade de ocorrer numa certa dura¸c˜ao (ou, de modo equivalente, at´e um determinado instante de tempo). Uma vez que o n´ucleo markoviano ´e a fun¸c˜ao respons´avel pela produ¸c˜ao de qualquer resultado deste modelo (da´ı a sua designa¸c˜ao), ´e comum utilizar PeFpara a constru¸c˜ao do n´ucleo Q, fazendo uso da rela¸c˜ao vista anteriormente. Com a finalidade de simplificar a constru¸c˜ao do modelo, considerem-se as fun¸c˜oes que definem: •A probabilidade de que o processo deixe um determinado estado iat´e ao instante t, sabendo que no instante stinha chegado ao estado i: Hi(s, t) = P(Tn+1 ≤t|Jn=i, Tn=s) = M X j=1 Qij(s, t). •a probabilidade de que o sistema chegue ao estado jno instante t, sabendo que no instante stinha chegado ao estado i: bij(s, t) = P(Jn+1 =j, Tn+1 =t|Jn=i, Tn=s). Note-se que esta probabilidade apenas pode ser definida quando se assume que o processo segue um percurso temporal discreto e que, por sua vez, pode ser escrita em fun¸c˜ao de Qij (s, t): bij(s, t) = Qij(s, t)−Qij(s, t −1) se s<t 0 se s=t Est˜ao agora definidas as ferramentas necess´arias para caracterizar um dos principais outputs deste modelo, e uma das suas fun¸c˜oes de interesse: o processo semi-markoviano de primeira ordem, n˜ao-homog´eneo em tempo discreto Z= (Z(t)), que representa, para cada instante de tempo t, o estado ocupado pelo processo. As probabilidades de transi¸c˜ao para este processo s˜ao dadas por φij(s, t) = P(Z(t) = j|Z(s) = i). Note-se que Z(t) e JN(t)s˜ao o mesmo processo. No entanto, φij(s, t) difere de Qij(s, t), na medida em que o primeiro considera a possibilidade de transi¸c˜oes interm´edias at´e `a chegada ao estado j, ao passo que o segundo considera a probabilidade de chegada a esse estado na transi¸c˜ao seguinte. 24 CAP´ ITULO 2. MODELAC¸ ˜ AO As probabilidades de transi¸c˜ao do processo Z(t) s˜ao definidas como φij(s, t) = δij(1 −Hi(s, t)) + M X β=1 t X ϑ=s biβ(s, ϑ)φβj(ϑ, t) [Janssen e Manca (2002)]. (2.1) onde δij representa o delta de Kronecker, dado por δij =1 se i=j 0 se i6=j.O resultado anterior foi provado pelos mesmos autores, que tamb´em demonstraram que esta equa¸c˜ao admite uma solu¸c˜ao ´unica [Janssen e Manca (2006)]. Para uma maior intui¸c˜ao do conceito, deve interpretar-se a express˜ao anterior: o primeiro termo representa a probabilidade de permanˆencia do processo no estado idurante os instantes set(porque o percurso do processo n˜ao ´e pr´e-definido e n˜ao existe garantia de que transite de estado) e o segundo termo representa a probabilidade de transi¸c˜ao do processo do estado ipara o estado jentre os instantes set(notando que entretanto pode passar por outros estados βem instantes interm´edios ϑ). Como esta ´e uma fun¸c˜ao diferencial que depende de Hi(s, t) e de biβ(s, ϑ), basta a defini¸c˜ao de uma condi¸c˜ao de fronteira para que se possa resolver. Quando s=t, e como no mesmo instante de tempo n˜ao se pode estar em dois estados diferentes, a condi¸c˜ao referida ser´a dada por φij(s, s) = 1 se i=j 0 se i6=j e a equa¸c˜ao pode ser resolvida do seguinte modo: •Define-se um valor m´aximo para o tempo que se vai estudar, Tmax. Define-se φij(s, s)∀s∈0, ..., Tmax; •Calcula-se φij(Tmax −1, Tmax). Analisando a equa¸c˜ao que define o processo, verifica-se que este ´e um c´alculo trivial, uma vez que s´o depende de φij(Tmax, Tmax), conhecido pelo passo anterior; •Calcula-se φij(Tmax −2, Tmax), que necessita de φij(Tmax −1, Tmax) e φij(Tmax, Tmax), ambos valores j´a obtidos; •Calcula-se φij(Tmax −3, Tmax), que depende de φij(Tmax −1, Tmax), φij(Tmax, Tmax) e φij(Tmax −2, Tmax), entretanto conhecidos; •Continua-se o processo at´e que φij(0, Tmax) esteja calculado; •Repete-se o processo, iniciando em φij(Tmax −1, Tmax−1). Uma vez resolvida esta equa¸c˜ao, ser´a poss´ıvel analisar todas as probabilidades de transi¸c˜ao e permanˆencia do sistema Z. Lembrando a defini¸c˜ao de N(t), seja ainda B(t) = t−TN(t),t∈N o tempo desde a ´ultima transi¸c˜ao, podendo ser interpretado como o processo temporal recursivo de Z(t) (backward process). Notando que (Z, B) ´e um processo de Markov, ´e agora 2.2. INTRODUC¸ ˜ AO AOS PROCESSOS SEMI-MARKOVIANOS 25 poss´ıvel escrever as diferentes fun¸c˜oes vistas anteriormente, mas tendo agora em considera¸c˜ao este novo processo B: bHi(u, s, t) = P(TN(s)+1 ≤t|JN(s)=i, TN(s)=u, TN(s)+1 > s) (u≤s<t) bQi,j(u, s, t) = P(TN(s)+1 ≤t, JN(s)+1 =j|JN(s)=i, TN(s)=u, TN(s)+1 > s) (u≤s<t) Observando que bHi(s, s, t) = Hi(s, t) e bQi,j(s, s, t) = Qi,j (s, t), conclui-se que bQi,j(u, s, t) = bQi,j(u, u, t) 1−bHi(u, u, t)=Qi,j(u, t) 1−Hi(u, s) Portanto bφi,j(u, s, t) = P(Z(t) = j|TN(s)=u, Z[u, s] = i) = =P(Z(t) = j|Z(s) = i, B(s) = s−u) devolve a probabilidade, sabendo que o sistema no instante de tempo sestava no estado ie que entrou nesse estado no instante u(por outras palavras, que estava no estado ih´a s−uunidades de tempo), de estar no instante tno estado j. Este ´e um resultado com muito potencial. Por exemplo, a probabilidade do sistema se encontrar num estado de vento considerado grave poder´a depender de h´a quanto tempo ´e que se encontra num estado com essa categoriza¸c˜ao. Para construir qualquer output do modelo ´e necess´ario estimar Q. H´a duas formas de o fazer: 1Estimar as fun¸c˜oes pij(s) e Fij(s, t). Foque-se a forma emp´ırica: definir pij(s) como a raz˜ao entre o n´umero de observa¸c˜oes que se encontravam no estado ino instante s cujo estado seguinte foi je o n´umero de observa¸c˜oes que se encontravam no estado ino instante s. Considerando aij como o n´umero de transi¸c˜oes que ocorreram do estado ipara o estado j,Barbu e Limnios (2009) mostram que o estimador emp´ırico ˆpij =aij PM j=1 aij se aproxima do estimador por m´axima verosimilhan¸ca. Ser´a este o utilizado nas simula¸c˜oes. Analogamente, definir empiricamente Fij(s, t) como a raz˜ao entre o n´umero de observa¸c˜oes que se encontravam no estado ino instante s, cujo estado seguinte foi je cuja transi¸c˜ao aconteceu no instante te o n´umero de observa¸c˜oes que se encontravam no estado ino instante s, cujo estado seguinte foi j. Observe-se ainda que pij(s) n˜ao tem necessariamente de depender de s. Pode assumir-se que pij(s) = pij(0) ∀s, pelo que a probabilidade de transi¸c˜ao entre estados seria a mesma ao longo do tempo. 2Estimar diretamente Qij(s, t). Para a estima¸c˜ao direta de Q, note-se que a estima¸c˜ao emp´ırica consistiria na an´alise da probabilidade de haver uma transi¸c˜ao no preciso instante de tempo em estudo, e fazer essa an´alise para todos os instantes de tempo considerados. Por ser menos intuitiva, apesar de exequ´ıvel, analisar-se-´a apenas o caso anterior. 32 CAP´ ITULO 3. PREVIS ˜ AO rede n˜ao ´e mais do que um conjunto de n´os, em que alguns est˜ao na camada de entrada (onde as unidades recebem os padr˜oes), alguns nas camadas interm´edias / escondidas (onde s˜ao efetuados o processamento e extra¸c˜ao de caracter´ısticas) e alguns na camada de sa´ıda (que apresenta o resultado final do processamento que se desencadeia na camada interm´edia). A l´ogica da rede neuronal est´a sobretudo nestes processamentos, que podem ser muito simples (cingindo-se `a soma de inputs, por exemplo) ou muito complexos (se um n´o contiver uma rede neuronal, por exemplo). Figura 3.2: Ilustra¸c˜ao de uma Rede Neuronal Artificial. Uma ANN tem a capacidade de aprender com a informa¸c˜ao que lhe ´e fornecida, melhorando o seu desempenho durante o processo de aprendizagem, uma vez que aprende por um processo iterativo de ajuste de pesos (for¸cas sin´apticas). Em cada itera¸c˜ao do processo de aprendizagem, apresenta a capacidade de aperfei¸coar a sua representa¸c˜ao porque aprende por treino (segundo certas regras pr´e-definidas). Ao conjunto de regras pr´e-definidas pelo qual se faz a altera¸c˜ao dos pesos d´a-se o nome de algoritmo de aprendizagem, que define a forma como os pesos s˜ao corrigidos e qual a estrutura da rede. O algoritmo mais utilizado ´e o Backpropagation. Por norma, um dos maiores problemas das ANN ´e a defini¸c˜ao da estrutura / arquitetura da rede, isto ´e, a estima¸c˜ao do n´umero de camadas ocultas e do n´umero de neur´onios em cada camada. Cada neur´onio recebe impulsos de entrada e calcula a informa¸c˜ao de sa´ıda como fun¸c˜ao desses impulsos, pelo que ´e realizado um c´alculo linear inicial nos inputs e, seguidamente, aplicada uma fun¸c˜ao de ativa¸c˜ao (Tabela 3.1). Fun¸c˜ao de Ativa¸c˜ao Express˜ao Linear f(x) = x Sinusoidal f(x) = sin(x) f(x) = cos(x) Sigm´oide f(x) = 2 1+e−x−1 Gaussiana f(x) = e−x2 2 Tabela 3.1: Fun¸c˜oes de Ativa¸c˜ao Usuais na Aplica¸c˜ao de Redes Neuronais Artificiais. Os sinais resultantes s˜ao posteriormente somados e `a soma resultante ´e aplicada uma fun¸c˜ao n˜ao linear - fun¸c˜ao transferˆencia - que verifica se o valor resultante da soma entre o produto 3.2. PREVIS ˜ AO - M´ ETODOS DE DATA MINING 33 dos sinais de entrada pelos respetivos pesos atingiu ou n˜ao um valor limite pr´e-definido, sendo assim gerado o output [Faraway (2005)]. Uma das arquiteturas de redes mais utilizada ´e a feedforward, tamb´em chamada rede sem realimenta¸c˜ao, que se caracteriza pelo agrupamento de neur´onios em camadas e pelo facto de o sinal percorrer a rede numa ´unica dire¸c˜ao (da entrada para a sa´ıda), n˜ao se estabelecendo liga¸c˜oes entre os neur´onios de uma mesma camada. Relembre-se ainda que, de acordo com um dos Teoremas da Aproxima¸c˜ao Universal [Suykens et al. (2012)], “uma rede feedforward com uma ´unica camada escondida que cont´em um n´umero finito de neur´onios, pode aproximar fun¸c˜oes cont´ınuas em subconjuntos compactos de Rn, sob suposi¸c˜oes leves na fun¸c˜ao de ativa¸c˜ao” [Gybenko (1989)]. Neste caso, a rede neuronal feedforward com uma camada escondida toma a forma fo(X h whofh(X i whixi)) onde forepresenta a fun¸c˜ao de transferˆencia, fhrepresenta a fun¸c˜ao de ativa¸c˜ao e w representa os pesos das liga¸c˜oes da rede (whonas liga¸c˜oes entre o input e os neur´onios da camada escondida e whinas liga¸c˜oes entre os neur´onios da camada escondida e o output). xirepresentam os inputs da rede. Ser´a esta a rede neuronal utilizada no caso pr´atico, a ver no cap´ıtulo de Aplica¸c˜ao. Observa¸c˜ao: Na estima¸c˜ao do n´umero de camadas e neur´onios da rede, dever˜ao ser tidos em conta os seguintes dois efeitos poss´ıveis: underfitting, caracterizado pela defini¸c˜ao de poucos neur´onios, n˜ao suficientes para que se consiga estabelecer uma rede neuronal fi´avel (ou para que sejam identificados padr˜oes) e overfitting, definida pela existˆencia de muitos neur´onios, que s˜ao treinados por um n´umero limitado de informa¸c˜ao contida no conjunto de dados. Uma rede neuronal treinada pode ser usada para fornecer proje¸c˜oes face a situa¸c˜oes de interesse. 3.2.2 ´ Arvores de Regress˜ao As ´arvores podem ser usadas em problemas de regress˜ao ou classifica¸c˜ao. A principal diferen¸ca reside no facto de as folhas das ´arvores de regress˜ao conterem previs˜oes num´ericas e n˜ao decis˜oes. O objetivo deste m´etodo consiste na parti¸c˜ao do espa¸co preditivo do conjunto de treino em regi˜oes, de modo que essas regi˜oes (subconjuntos finais) sejam t˜ao “puras” quanto poss´ıvel. A parti¸c˜ao passo a passo obtida corresponde a uma aproxima¸c˜ao da parti¸c˜ao ´otima. Para cada n´o da ´arvore, ´e necess´ario escolher a vari´avel que melhor segmenta esse n´o, definindo-se uma medida de impureza de tal modo que os descendentes do n´o sejam mais puros (existindo menos mistura de informa¸c˜ao das regi˜oes at´e a´ı definidas) do que o n´o que lhes deu origem. 34 CAP´ ITULO 3. PREVIS ˜ AO No caso de uma ´arvore de classifica¸c˜ao, as medidas de impureza para um dado n´o que s˜ao mais frequentemente utilizadas (por favorecem os n´os mais puros quando comparadas com o erro de classifica¸c˜ao), s˜ao: •Quantidade de Informa¸c˜ao de Shannon: −PjP(cj)log2(P(cj)) •´ Indice de Gini: 1 −PjP2(cj) Note-se que P(cj) ´e a fra¸c˜ao das vari´aveis independentes no n´o em an´alise que pertencem `a regi˜ao cj(sejam c= (c1, c2, ...) as regi˜oes definidas pelo algoritmo de aprendizagem). Estas medidas satisfazem as propriedades definidas em Breiman et al. (1984). No caso de uma ´arvore de regress˜ao, o custo de escolher o valor y=anum dado n´o ´e em geral determinado por uma das duas seguintes medidas: E[(Y−a)2] ou E[|Y−a|]. A a¸c˜ao que minimiza o custo, no primeiro caso, ´e a atribui¸c˜ao a ado valor da m´edia de Y(a=E[Y]), enquanto que no segundo caso ´e a atribui¸c˜ao a ado valor da mediana de Y (a=εY ). Por este motivo, as medidas de impureza a considerar em regress˜ao s˜ao: •Desvio quadr´atico m´edio: E[(Y−E[Y])2] •Desvio absoluto m´edio: E[|Y−εY |] A medida de impureza utilizada no caso pr´atico ser´a a dada pelo desvio quadr´atico m´edio. 3.2.3 Regress˜ao com M´aquinas de Suporte Vetorial A no¸c˜ao de Regress˜ao com M´aquinas de Suporte Vectorial (Support Vector Regression, SVRs) [Vapnik (1995)] surgiu como uma generaliza¸c˜ao das M´aquinas de Suporte Vetorial (SVMs), que por sua vez surgiram para resolver os problemas computacionais do m´etodo do n´ucleo. Estes estavam relacionados com a lentid˜ao de resposta para grandes conjuntos de dados e necessidade de armazenamento de toda a base de dados para fazer previs˜ao. De modo informal, o m´etodo do n´ucleo ´e um m´etodo n˜ao param´etrico para estima¸c˜ao de curvas de densidades onde cada observa¸c˜ao ´e ponderada pela distˆancia em rela¸c˜ao a um valor central, designado como n´ucleo, cuja ideia foi introduzida para previs˜ao por Nadaraya (1964) e Watson (1964). O m´etodo SVR ´e complexo e a sua explica¸c˜ao n˜ao ser´a exaustiva, mas intuitiva (ver Figura 3.3). A explora¸c˜ao deste m´etodo pode ser consultada, por exemplo, no tutorial Burges (1998). De um modo muito geral, nas SVMs, atrav´es da aplica¸c˜ao de uma fun¸c˜ao de n´ucleo (K) `as observa¸c˜oes, estas s˜ao projetadas num espa¸co de maior dimens˜ao no qual os dados podem ser separados por um hiperplano. Quando os dados de treino s˜ao separ´aveis, o hiperplano ´otimo no espa¸co caracter´ıstico apresenta a m´axima margem de separa¸c˜ao. 3.3. AVALIAC¸ ˜ AO DOS MODELOS 35 Figura 3.3: Ilustra¸c˜ao de uma SVM. Dado um conjunto de treino (xi, yi)i, onde xi∈Rney∈ {1,−1}l, procura-se resolver o seguinte problema de otimiza¸c˜ao: min w, b, ξ 1 2wTw+C l X i=1 ξi sujeito a yi(wTκ(xi) + b)≥1−ξii= 1, . . . , l. ξi≥0, onde xirepresenta os vetores de treino, wTrepresenta a transposta de w,C > 0 o parˆametro de penaliza¸c˜ao dos erros ξieK(xi,xj)≡κ(xi)Tκ(xj) ´e a chamada fun¸c˜ao n´ucleo (ver Tabela 3.2). A fun¸c˜ao n´ucleo utilizada ser´a a radial. Fun¸c˜ao N´ucleo Express˜ao Linear xi Txj Polinomial (γxi Txj+r)dγ > 0 Radial exp(−γ||xi Txj||2)γ > 0 Sigm´oide tanh(γxi Txj+r) Tabela 3.2: Fun¸c˜oes N´ucleo Usuais na Aplica¸c˜ao de SVMs. 3.3 Avalia¸c˜ao dos Modelos Para fazer a avalia¸c˜ao dos modelos de previs˜ao, o conjunto de dados ´e particionado em dois subconjuntos: um de treino, usado para a aprendizagem do modelo, e outro de teste, para analisar o erro associado ao modelo. Se assim n˜ao fosse, o modelo encontrado j´a estaria ajustado aos dados de teste, produzindo previs˜oes otimistas (e falseando, de algum modo, os resultados). A separa¸c˜ao pode ser realizada por valida¸c˜ao cruzada (m´etodo k-fold [Witten e Frank (2005)]) no caso dos algoritmos que aplicam os m´etodos de Data Mining acima referidos. Por´em, ´e muito importante lembrar que, para s´eries temporais, qualquer forma de reamostragem altera a ordem inicial dos dados e, por isso mesmo, esta t´ecnica n˜ao pode ser aplicada, pelo que, nesse caso em particular, foi usado um m´etodo designado por Sliding Window, que consiste na divis˜ao dos dados existentes em duas janelas: uma que cont´em as observa¸c˜oes 36 CAP´ ITULO 3. PREVIS ˜ AO anteriores a um dado instante de tempo e outra que cont´em as restantes observa¸c˜oes, e onde, para a aprendizagem no conjunto de teste, ´e constru´ıdo um novo modelo, para cada conjunto de teste, obtido treinando todos os dados anteriores a ele (note-se que, de cada vez que uma nova observa¸c˜ao ´e adicionada ao conjunto de treino, uma mais antiga ´e removida). As principais m´etricas usadas para avaliar os erros de previs˜ao dos modelos anteriormente explicados s˜ao definidas na Tabela 3.3, onde nrepresenta o n´umero de observa¸c˜oes dispon´ıveis, Yto valor observado em teˆ Yta previs˜ao do modelo para esse mesmo valor. Denomina¸c˜ao Acr´onimo Express˜ao Erro Quadr´atico M´edio MSE 1 nPn i=1(Yt−ˆ Yt)2 Ra´ız do Erro Quadr´atico M´edio RMSE √MSE Erro Absoluto M´edio MAD 1 nPn i=1 |Yt−ˆ Yt| Erro Percentual Absoluto M´edio MAPE 100 nPn i=1 |Yt−ˆ Yt Yt| Tabela 3.3: M´etricas Comuns na Avalia¸c˜ao dos Erros dos Modelos de Previs˜ao. Observe-se que o MSE apresenta uma grande sensibilidade a erros elevados, uma vez que considera o quadrado das diferen¸cas entre os valores observados e previstos. ´ E por esse motivo que muitas vezes se considera o RMSE, que atenua essa desvantagem. O MAD, quando comparado com o MSE, apresenta a vantagem de avaliar os erros na unidade original dos dados (e n˜ao ao quadrado), tratando-os de igual modo. O MAPE, que ser´a o erro usado no caso pr´atico, mede o “tamanho do erro” e corresponde `a m´etrica que mais informa¸c˜ao d´a relativamente `a qualidade preditiva do modelo. Note-se que, nesta ´ultima m´etrica, as observa¸c˜oes tais que Yt= 0 n˜ao podem ser avaliadas. Uma vez definidas as ferramentas te´oricas fundamentais propostas para utiliza¸c˜ao, ´e poss´ıvel continuar a an´alise, agora sobre a aplica¸c˜ao das mesmas a dados reais. FCUP 37 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Cap´ıtulo 4 Aplica¸c˜ao A base de dados analisada apresenta algumas vari´aveis clim´aticas (referidas na Tabela 4.1 e cuja nomenclatura ser´a, daqui em diante, a presente na segunda coluna da mesma) cujas caracter´ısticas s˜ao detalhadas no Anexo D. Vari´avel Abreviatura Unidade de Medida Ano A - Mˆes M - Dia D - Hora H - Esta¸c˜ao do Ano E - Dire¸c˜ao do Vento Hor´aria DV Graus Precipita¸c˜ao Total Hor´aria PT mm Temperatura M´ınima Hor´aria TMin Graus Celsius Temperatura M´edia Hor´aria TMed Graus Celsius Temperatura M´axima Hor´aria TMax Graus Celsius Humidade Relativa M´ınima Hor´aria HMin % Humidade Relativa M´edia Hor´aria HMed % Humidade Relativa M´axima Hor´aria HMax % Radia¸c˜ao Global M´ınima Hor´aria RMin W/m2 Radia¸c˜ao Global M´edia Hor´aria RMed W/m2 Radia¸c˜ao Global M´axima Hor´aria RMax W/m2 Velocidade de Vento M´ınima Hor´aria VMin m/s Velocidade de Vento M´edia Hor´aria VMed m/s Velocidade de Vento M´axima Hor´aria VMax m/s Tabela 4.1: Vari´aveis da Base de Dados. Para se avaliar quais as vari´aveis com maior influˆencia sobre a velocidade m´axima de vento, foi efetuada uma pequena an´alise, precedida de um tratamento da base de dados. Neste ´ultimo, executou-se : •preenchimento de valores em falha (vulgarmente designados por NA - Not Available), que surgiam em seis das vari´aveis (HMed, HMin e HMax, em aproximadamente 3% 38 CAP´ ITULO 4. APLICAC¸ ˜ AO do tamanho da amostra; PT, em duas observa¸c˜oes; VMin e VMax, numa observa¸c˜ao), substituindo-os pelo correspondente valor da mediana da vari´avel (seguindo a abordagem sugerida por Torgo (2010)). •dada a eventual importˆancia que a vari´avel temporal “Data” poder´a ter no estudo, esta foi substitu´ıda por cinco outras vari´aveis: Ano, Mˆes, Dia, Hora e Esta¸c˜ao do Ano. Tal como escrito anteriormente, o foco incidir´a sobre a vari´avel vento, em particular sobre a vari´avel velocidade m´axima do vento da base de dados. Como as instru¸c˜oes implementadas em R para o c´alculo da correla¸c˜ao n˜ao aceitam vari´aveis categ´oricas, a vari´avel Esta¸c˜ao do Ano foi tranformada em num´erica apenas para que este c´alculo pudesse ter em conta todas as vari´aveis da base de dados. Note-se que, dado um conjunto de dados com diversas vari´aveis, se o objetivo for diminuir o seu n´umero, n˜ao se deve considerar as que est˜ao altamente correlacionadas entre si (uma vez que contˆem a mesma informa¸c˜ao, tornando-se redundantes). Pela observa¸c˜ao de boxplots e pela matriz de correla¸c˜ao entre vari´aveis, n˜ao aparenta existir uma dependˆencia clara entre a vari´avel em estudo e as restantes (exceto no caso em que a vari´avel clim´atica ´e a mesma, mas analisada em rela¸c˜ao aos seus valores m´aximo, m´ınimo e m´edio hor´arios). Figura 4.1: Correla¸c˜oes Amostrais. No entanto, o facto de n˜ao aparentar existir correla¸c˜ao linear entre as vari´aveis, quando analisadas duas a duas, n˜ao significa que a vari´avel em estudo n˜ao possa ser descrita por uma combina¸c˜ao das outras vari´aveis. Nesse caso, poder-se-ia proceder ao estudo dessa rela¸c˜ao considerando a velocidade m´axima de vento como vari´avel resposta de um modelo de regress˜ao linear m´ultipla e essa vari´avel, seguidamente designada por V, poderia ser escrita como V=β0+Xβ1+... +βpXp+u=Xβ onde X1, ..., Xps˜ao as vari´aveis explicativas, uos erros (tamb´em designados por res´ıduos) do modelo e β= (β0, β1, ...., βp) os coeficientes da regress˜ao a estimar. 39 No modelo cl´assico de regress˜ao linear assume-se que os erros us˜ao i.i.d., u∼N(0; ε=σ2Id) eVsegue uma distribu¸c˜ao normal com m´edia dependente de modo linear das vari´aveis explicativas, isto ´e, V|X∼N(µ(X), σ2(X)), E(V|X) = β0+Xβ1+... +βpXp. De um modo sucinto, face `a express˜ao geral de defini¸c˜ao do modelo de regress˜ao linear e para que se fa¸ca a correta interpreta¸c˜ao dos coeficientes estimados no modelo, ´e importante referir que: •O coeficiente independente, β0, representa a resposta caso todas as vari´aveis explicativas sejam nulas. •Para j= 1, ..., p o coeficiente jrepresenta o incremento m´edio de Vquando a vari´avel explicativa Xj´e aumentada de uma unidade e as restantes vari´aveis explicativas se mantˆem constantes, o que permite avaliar a intensidade da rela¸c˜ao entre VeXj. A primeira avalia¸c˜ao dos modelos deve ser feita por an´alise de testes de hip´oteses, crit´erios de informa¸c˜ao eres´ıduos. Antes de se observarem os resultados, ´e necess´ario salientar alguns aspetos: •Como o problema a resolver ´e o da estima¸c˜ao dos parˆametros βeσ2, mantendo σ2 fixo, βpode ser estimado pelo m´etodo da m´axima verosimilhan¸ca, pelo que, supondo a independˆencia entre as observa¸c˜oes, a fun¸c˜ao de verosimilhan¸ca a maximizar ser´a dada por L(β, σ2|(vi, xi)i) = n Y i=1 1 √2πσ2exp−1 2 (vi−ˆvi)2 σ2. Como a fun¸c˜ao logaritmo ´e crescente, maximizar Lequivale a maximizar log(L), ou seja l(β, σ2|(vi, xi)i) = log(L(β, σ2|(vi, xi)i)) =−n 2log(2πσ2)−1 2 n X i=1 (vi−ˆvi)2 σ2=−n 2log(2πσ2)−1 2σ2 n X i=1 (vi−Xi.β)2. Os crit´erios de informa¸c˜ao s˜ao usados para compara¸c˜ao de modelos n˜ao encaixados e aplicados a modelos constru´ıdos a partir da maximiza¸c˜ao do logaritmo da verosimilhan¸ca acima definido. Caracterizam-se pela penaliza¸c˜ao de modelos com maior n´umero de parˆametros e os mais usados s˜ao, segundo Pinheiro e Bates (2000): –crit´erio de informa¸c˜ao de Akaike (AIC), dado por −2l(ˆ β, ˆσ2)+2npar –crit´erio de informa¸c˜ao Bayesiana (BIC), dado por −2l(ˆ β, ˆσ2)+2nparlog(N) onde npar ´e o n´umero de parˆametros do modelo e No n´umero de observa¸c˜oes. 40 CAP´ ITULO 4. APLICAC¸ ˜ AO •Os principais testes envolvidos no teste de um modelo de regress˜ao linear s˜ao os seguintes [Faraway (2009)]: –teste de Wald, que indica a significˆancia para o modelo de cada um dos coeficientes estimados, testando a hip´otese nula de que βj= 0 para algum j= 0,1, ..., p; a an´alise dos resultados deste teste deve ser feita com precau¸c˜ao, uma vez que a rejei¸c˜ao da hip´otese nula n˜ao implica que todas as vari´aveis n˜ao significativas devam ser exclu´ıdas do modelo. A remo¸c˜ao de uma vari´avel explicativa pode fazer com que outras, que antes n˜ao eram consideradas significativas, o passem a ser. –teste-F, atrav´es do qual se pretende averiguar se um determinado grupo de vari´aveis explicativas ´e ou n˜ao significativo na explica¸c˜ao da vari´avel resposta do modelo de regress˜ao, testando a hip´otese nula de que β0=β1=... = βp= 0, e a alternativa dada por ∃j∈ {1, ..., p}tal que βj6= 0. ´ E usual apresentar o resultado de uma an´alise de regress˜ao sob a forma de uma tabela de an´alise da variˆancia, ANOVA, onde se indicam alguns valores necess´arios ao desenvolvimento do teste de hip´oteses anterior. Al´em disso, a varia¸c˜ao da vari´avel de resposta ´e decomposta na soma da varia¸c˜ao devida `a regress˜ao (Regression Sum of Squares) e da varia¸c˜ao residual (Residual Sum of Squares). Por esse motivo, tal como intuitivamente faz sentido, um modelo que apresente um bom ajustamento aos dados ser´a um modelo em que a varia¸c˜ao total ser´a essencialmente devida `a regress˜ao, apresentando uma varia¸c˜ao residual baixa. Neste ˆambito, foram avaliados modelos com e sem algumas das vari´aveis. As referentes aos valores m´edios, uma vez que correspondem a valores calculados e n˜ao medidos, n˜ao foram consideradas. Como visto na Figura 4.1, a informa¸c˜ao recolhida mostrou que algumas das vari´aveis s˜ao altamente correlacionadas. Para decidir quais dessas vari´aveis se deveria incluir no modelo, optou-se pelas menos correlacionadas entre si. Consequentemente, retiveramse as vari´aveis TMin, HMax e RMin em detrimento das vari´aveis TMax, HMin e RMax, respetivamente. As diferentes aproxima¸c˜oes do modelo por regress˜ao linear sugeriram que algumas das vari´aveis n˜ao eram significativas (valor-p superior a 0.05). Estas foram removidas uma a uma e os resultados analisados. O resultado final dessa an´alise, obtido em R, foi o seguinte: Call: lm(formula = VMax ~ I(A) + I(M) + I(Estacao) + TMin + HMax + RMin + DV + PT, data = BaseDados) Residuals: Min 1Q Median 3Q Max -16.0302 -1.9165 -0.2585 1.6601 18.3484 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -1.941e+02 9.442e+01 -2.056 0.0398 * I(A) 1.025e-01 4.686e-02 2.188 0.0287 * I(M) 2.951e-02 1.328e-02 2.222 0.0263 * I(Estacao)Outono 4.924e-01 1.242e-01 3.966 7.35e-05 *** 41 I(Estacao)Primavera -8.014e-01 8.810e-02 -9.097 < 2e-16 *** I(Estacao)Ver~ao -3.265e-01 1.288e-01 -2.535 0.0112 * TMin -1.795e-01 7.098e-03 -25.294 < 2e-16 *** HMax -3.998e-02 1.519e-03 -26.330 < 2e-16 *** RMin 2.204e-03 1.537e-04 14.341 < 2e-16 *** DV -1.347e-02 2.179e-04 -61.822 < 2e-16 *** PT 4.609e-01 2.431e-02 18.955 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.227 on 17732 degrees of freedom Multiple R-squared: 0.3568,Adjusted R-squared: 0.3564 F-statistic: 983.4 on 10 and 17732 DF, p-value: < 2.2e-16 Coment´ario: A rejei¸c˜ao da hip´otese nula no teste-F n˜ao significa que o modelo seja bom. Ali´as, observando o valor do coeficiente de determina¸c˜ao (quadrado do coeficiente de correla¸c˜ao linear de Pearson amostral entre os valores observados e os ajustados), que deveria ser pr´oximo de 1, pode concluir-se que o modelo apresenta uma varia¸c˜ao residual alta (o coeficiente de determina¸c˜ao ´e igual a 0.3568), pelo que ´e um indicador de um mau ajustamento aos dados. A inexistˆencia de rela¸c˜ao linear entre a vari´avel resposta e as vari´aveis explicativas conduz a valores de coeficiente de determina¸c˜ao pr´oximos de 0. Este tamb´em n˜ao parece ser o caso. Vai ser necess´ario analisar os res´ıduos e outras medidas de significˆancia. Figura 4.2: An´alise Gr´afica de Res´ıduos (caso hor´ario). Chegou-se `a conclus˜ao que o modelo apresenta 2485 observa¸c˜oes com res´ıduo superior a 3.3 que, segundo a literatura da ´area, ´e o valor de referˆencia a partir do qual as observa¸c˜oes devem ser analisadas uma a uma, com o intuito de perceber qual a sua influˆencia sobre o modelo estimado pela regress˜ao, o que n˜ao parece ser exequ´ıvel dado o elevado n´umero de observa¸c˜oes nessas condi¸c˜oes! 48 CAP´ ITULO 4. APLICAC¸ ˜ AO Com base em todo o processo explicado na sec¸c˜ao te´orica de modela¸c˜ao desta tese, foi obtida a matriz de transi¸c˜oes do sistema (ver equa¸c˜ao (2.1) referida na Introdu¸c˜ao aos Processos Semi-Markovianos). A explora¸c˜ao da mesma poderia ser num´erica, mas optou-se pela an´alise gr´afica, uma vez que permite uma leitura mais intuitiva acerca do potencial da informa¸c˜ao obtida por este modelo. A vari´avel estudada foi discretizada em 6 estados: dois de vento fraco, dois de vento interm´edio e dois de vento forte (ver Tabela 3do Anexo C). Alguns dos gr´aficos obtidos pelo estudo de φij(1, t), t = 1, ..., Tmax (efetuado a partir da primeira hora registada na base de dados, para Tmax = 50 horas) foram os seguintes: Para i= 1 Para i= 2 Para i= 3 Para i= 4 Para i= 5 Para i= 6 Tabela 4.4: Probabilidades de Transi¸c˜ao Estimadas (da Cadeia Semi-markoviana de Primeira Ordem). Coment´ario: A an´alise gr´afica permite constatar que, no final das 50 horas analisadas, existe uma estabiliza¸c˜ao dos valores das probabilidades de transi¸c˜ao. ´ E tamb´em poss´ıvel observar que, nos estados de vento fraco (1 e 2), a probabilidade de transi¸c˜ao para estado de vento forte s˜ao muito baixas e que o acontecimento mais prov´avel ´e a transi¸c˜ao para o 4.1. MODELAC¸ ˜ AO 49 outro estado de vento fraco. Quanto aos estados de vento forte (5 e 6), ´e interessante verificar que, ap´os algumas horas, ´e mais prov´avel que o sistema retorne a um estado de vento fraco do que se mantenha num de vento forte ou interm´edio (3 e 4). Este facto permite, de algum modo, confirmar o que o senso comum nos diz: “depois da tempestade vem a bonan¸ca”. No caso do estado 5, nas primeiras horas, ´e mais prov´avel que transite para um estado interm´edio (3 ou 4), ao passo que no estado 6 ´e mais prov´avel que transite para um estado de vento forte (5). As primeiras 15 horas tamb´em parecem ser determinantes no decorrer do processo, em particular nos estados de vento interm´edios e nos estados de vento forte, na medida em que a probabilidade de transi¸c˜ao para estados com velocidade de vento pr´oximas das que o estado original apresenta s˜ao mais prov´aveis nessas primeiras horas do que nas seguintes. Note-se que s´o se observou o comportamento do sistema nas primeiras horas uma vez que a matriz de transi¸c˜ao Pconsiderada ´e igual para todos os instantes de tempo. Se fosse diferente para cada instante, o estudo de outros dom´ınios temporais poderia ser mais informativo. Esta fun¸c˜ao fornece quase toda a informa¸c˜ao necess´aria sobre o processo, uma vez que ´e poss´ıvel ao leitor estudar as probabilidades de transi¸c˜ao entre os estados de interesse, nos instantes de interesse. A explora¸c˜ao desta informa¸c˜ao seria ´util, por exemplo, caso existissem v´arias regi˜oes em an´alise. Nesse caso, seria poss´ıvel a uma qualquer companhia seguradora analisar qual o comportamento esperado do vento nas diferentes regi˜oes, fazendo diferentes atribui¸c˜oes de pr´emios (pricing) consoante a regi˜ao do segurado. Do mesmo modo, esta informa¸c˜ao poderia ser usada para estudar o potencial energ´etico de uma determinada regi˜ao (acrescentando ao estudo outro tipo de vari´aveis), aliando informa¸c˜ao relativa n˜ao s´o `a velocidade de vento, mas tamb´em `a sua dire¸c˜ao [D’Amico et al. (2012)]. Fez-se uma an´alise de diagn´ostico com base num histograma para se perceber se este seria semelhante `a s´erie original de forma a testar empiricamente se o c´alculo da matriz Ffoi feito corretamente e se as cadeias semi-markovianas (de primeira e segunda ordens) conseguiam modelar os dados ou n˜ao. Os resultados obtidos relativos `as simula¸c˜oes foram os seguintes: Figura 4.5: Simula¸c˜oes e respetivas Fun¸c˜oes de Correla¸c˜ao. Coment´ario: Em qualquer um dos gr´aficos, os processos semi-markovianos aparentam fazer uma modela¸c˜ao que se aproxima do processo original de um modo mais coerente do que o processo markoviano. Note-se que o gr´afico da esquerda apenas permite avaliar a 50 CAP´ ITULO 4. APLICAC¸ ˜ AO representatividade dos estados no final do processo, mostrando o n´umero de vezes que o sistema esteve em cada um dos estados. A observa¸c˜ao do gr´afico de autocorrela¸c˜ao (lado direito da Figura 4.5) tamb´em corrobora a primeira afirma¸c˜ao, na medida em que a an´alise temporal definida por esta medida sugere um comportamento por parte dos processos semimarkovianos mais pr´oximo do original, quando comparada com a autocorrela¸c˜ao do processo markoviano. No entanto, a an´alise gr´afica n˜ao basta para fazer conclus˜oes. Para resultados mais elucidativos, mostram-se tabelas com quantidades de interesse, resultantes de 1000 simula¸c˜oes Monte Carlo dos trˆes processos. A Tabela 4.5 exibe as probabilidades m´edias de transi¸c˜ao estimadas (medidas em percentagem e arredondadas a duas casas decimais) e a Tabela 4.6 os tempos m´edios de espera estimados (hor´arios, arredondados a duas casas decimais) antes de ocorrer transi¸c˜ao. Os resultados mais pr´oximos dos originais (avaliados pela diferen¸ca em valor absoluto) est˜ao assinalados a azul 3. Reais Markov 123456 1−96.28 3.72 0.00 0.00 0.00 2 64.12 −34.67 1.21 0.00 0.00 3 1.34 60.22 −37.10 1.12 0.22 4 0.23 4.22 77.00 −18.08 0.47 5 0.00 1.14 14.77 80.68 −3.41 6 0.00 0.00 14.28 71.43 14.29 − 123456 1−95.64 4.36 0.00 0.00 0.00 262.46 −35.89 1.65 0.00 0.00 3 1.91 63.74 −33.06 1.29 0.00 4 0.16 6.13 78.39 −14.84 0.48 5 0.00 0.00 16.97 78.57 −4.46 6 0.00 0.00 12.50 75.00 12.50 − Semi-Markov de Ordem 1 Semi-Markov de Ordem 2 1 2 3 4 5 6 1−95.86 4.14 0.00 0.00 0.00 2 65.59 −33.10 1.31 0.00 0.00 31.24 57.80 −39.28 1.57 0.11 4 0.44 2.44 78.93 −17.52 0.67 5 0.00 1.08 15.05 79.57 −4.30 6 0.00 0.00 0.00 100.00 0.00 − 123456 1−96.40 3.60 0.00 0.00 0.00 2 65.52 −33.44 1.04 0.00 0.00 3 1.48 60.27 −36.77 1.25 0.23 40.25 5.42 78.08 −15.76 0.49 5 0.00 2.59 12.99 84.42 −0.00 6 0.00 0.00 25.00 25.00 50.00 − Tabela 4.5: Probabilidades M´edias de Transi¸c˜ao Estimadas (em %). Coment´ario: A informa¸c˜ao presente na primeira tabela n˜ao indica que exista efetivamente um melhor ajuste dos modelos semi-markovianos. Apesar de grande parte dos melhores resultados estar assinalada sobre os resultados da modela¸c˜ao semi-markoviana, ainda h´a situa¸c˜oes em que o modelo markoviano ´e o melhor. De facto, existem situa¸c˜oes em que os processos semi-makovianos se aproximam dos originais (como por exemplo nas transi¸c˜oes 1→2, 2 →4, 4 →3 e 5 →6) em detrimento do de Markov. H´a outras situa¸c˜oes em que esse facto n˜ao se verifica (como por exemplo nas transi¸c˜oes do estado 6 para qualquer um dos outros estados). Ser´a necess´ario um crit´erio de decis˜ao mais forte para decidir qual o processo que melhor aproxima as caracter´ısticas do processo original. 3caso existam valores iguais, s˜ao todos assinalados. Para esta avalia¸c˜ao s˜ao consideradas apenas as probabilidades e tempos diferentes de 0. 4.1. MODELAC¸ ˜ AO 51 Por esse motivo, foram analisados os tempos de espera, onde os modelos semi-markovianos apresentam claramente melhores resultados do que o modelo markoviano, em particular nos estados de vento fraco e interm´edio. Reais Markov 123456 1−8.56 9.44 0.00 0.00 0.00 2 3.42 −3.10 4.06 0.00 0.00 3 1.33 3.03 −3.21 1.60 4.50 4 1.00 1.28 2.12 −2.56 1.00 5 0.00 1.00 1.08 1.80 −1.00 6 0.00 0.00 1.00 1.00 1.00 − 123456 1−3.66 3.65 0.00 0.00 0.00 2 2.83 −2.83 2.80 0.00 0.00 3 2.44 2.45 −2.45 2.43 0.00 4 1.82 1.81 1.82 −1.83 1.83 5 0.00 1.30 1.30 1.29 −1.30 6 0.00 0.00 1.00 1.00 1.00 − Semi-Markov de Ordem 1 Semi-Markov de Ordem 2 123456 1−8.55 9.47 0.00 0.00 0.00 2 3.43 −3.10 4.07 0.00 0.00 31.33 3.03 −3.20 1.60 4.42 41.00 1.28 2.12 −2.57 1.00 5 0.00 1.00 1.08 1.81 −1.00 6 0.00 0.00 1.00 1.00 1.00 − 123456 1−8.55 9.43 2.00 0.00 0.00 23.42 −3.10 4.03 0.00 0.00 3 1.34 3.03 −3.21 1.60 4.56 41.00 1.28 2.12 −2.56 1.00 5 0.00 1.00 1.07 1.80 −1.00 6 0.00 0.00 1.00 1.00 1.00 − Tabela 4.6: Tempos de Espera M´edios Estimados. Esta an´alise tamb´em permite inferir alguma informa¸c˜ao interessante sobre o processo de velocidade de ventos na regi˜ao de Settala: as transi¸c˜oes entre estados de vento forte ocorrem num intervalo temporal curto, uma vez que o tempo m´edio de espera para a ocorrˆencia de transi¸c˜ao do estado 5 para o estado 6 ou do estado 6 para o estado 5 ´e de uma hora. Em m´edia, o tempo que o sistema demora a transitar para o estado 2, quando vem do estado 1, ´e de aproximadamente 8.56 horas. Curiosamente, ´e mais r´apido a transitar para o estado 1 se provier do estado 2 - demora aproximadamente 3.43 horas. No entanto, relembrando os crit´erios climatol´ogicos usados para a discretiza¸c˜ao em estados, ´e natural que os tempos de transi¸c˜ao nos estados de vento fraco tenham uma banda temporal mais larga, uma vez que contˆem velocidades de vento baixas e comuns na regi˜ao, que podem variar de modo gradual ao longo do dia (note-se que varia¸c˜oes de 1 m/s para 3 m/s ou de 4.5 m/s para 8.5 m/s, por exemplo, est˜ao contempladas no mesmo estado - estados 1 ou 2, respetivamente - e por isso n˜ao s˜ao controladas). Em rela¸c˜ao aos resultados temporais, realce-se que eram expect´aveis, uma vez que a grande diferen¸ca entre a abordagem tradicional (de Markov) e a realizada na presente tese ´e a da considera¸c˜ao de informa¸c˜ao temporal. Antes de terminar a an´alise da modela¸c˜ao, deve salientar-se que, em rela¸c˜ao `as quantidades de interesse avaliadas nesta sec¸c˜ao, n˜ao foram observadas melhorias significativas pela considera¸c˜ao do processo semi-markoviano de segunda ordem em rela¸c˜ao ao de primeira ordem. No entanto, como j´a foi referido, o processo considerado ´e de segunda ordem nos estados e de primeira ordem no tempo. D’Amico et al. (2013) referem melhorias na 52 CAP´ ITULO 4. APLICAC¸ ˜ AO aproxima¸c˜ao a dados reais por modelos semi-markovianos de segunda ordem nos estados e no tempo. Globalmente, os resultados vistos permitem afirmar que a modela¸c˜ao semi-markoviana aparenta fazer um melhor ajuste aos dados reais de velocidade m´axima de vento hor´aria do que a modela¸c˜ao markoviana. 4.2 Previs˜ao - Mecanismos de Explora¸c˜ao Relativamente aos processos de previs˜ao de velocidade de vento, face aos maus resultados obtidos nas tentativas de previs˜ao hor´aria, a partir deste momento a an´alise ser´a feita sobre os dados di´arios m´aximos de velocidade de vento. A s´erie a analisar ´e a da Figura 4.6. Figura 4.6: S´erie Temporal Di´aria. Figura 4.7: ARIMA(3,0,1). Tentou modelar-se a s´erie temporal por um ARIMA usando a instru¸c˜ao auto.arima do R (cujo algoritmo foi desenvolvido, apresentado e implementado em Hyndman e Khandakar (2007)). Esta fun¸c˜ao sugere o modelo ARIMA que melhor se adapta aos dados (com base em testes de ra´ız unit´aria para determinar o valor de de com base nos crit´erios AIC, BIC ou AICc (AIC corrigido) para os outros dois parˆametros). Segundo os testes de Dickey-Fuller e Phillips-Perron referidos na parte te´orica, ´e rejeitada a hip´otese de existir pelo menos uma ra´ız dentro do c´ırculo unit´ario, para um n´ıvel de significˆancia de 5%. O teste de KPSS n˜ao rejeita a hip´otese de estacionariedade da s´erie. Os valores-p para cada teste foram, respetivamente, 0.01, 0.01 e 0.1. No caso destes testes, n˜ao ser´a poss´ıvel haver concordˆancia entre os resultados, uma vez que a hip´otese de interesse ´e a hip´otese nula num dos testes e a alternativa no outro. Como a n˜ao rejei¸c˜ao de uma hip´otese nula n˜ao implica a sua aceita¸c˜ao, pode dizer-se apenas que, segundo os dois primeiros testes, a hip´otese de n˜ao estacionariedade ´e rejeitada e segundo o terceiro teste, a estacionariedade n˜ao ´e rejeitada. Deste modo, sup˜oe-se que o modelo ARIMA poder´a ser aplicado a esta s´erie sem necessidade de a diferenciar. 4.2. PREVIS ˜ AO - MECANISMOS DE EXPLORAC¸ ˜ AO 53 Face `a previs˜ao por s´eries temporais, note-se que a reta de referˆencia no correlograma da ACF ´e dada por y=−1 n±(2 √n), e portanto depende do n´umero de observa¸c˜oes (n). Em seguida apresentam-se os correlogramas das autocorrela¸c˜oes totais e parciais, respetivamente, da s´erie temporal em an´alise. Figura 4.8: ACF. Figura 4.9: ACF (500 lags). Figura 4.10: PACF. A ACF n˜ao estabiliza abaixo da reta de referˆencia nos primeiros lags observados. Para um maior n´umero de lags, a ACF parece ter um comportamento an´alogo ao de uma onda criticamente amortecida. Como nem sempre a an´alise gr´afica ´e correta e suficiente para justificar uma decis˜ao, aliada ao facto de n˜ao se encontrarem referˆencias com um suporte te´orico s´olido aos m´etodos gr´aficos de apoio `a decis˜ao do modelo, devem testar-se v´arias ordens e escolher a que apresentar melhores resultados. Neste caso, como os dois modelos referidos nesta sec¸c˜ao foram utilizados sobretudo como mecanismos de explora¸c˜ao da s´erie temporal, esse estudo n˜ao foi feito. O modelo escolhido pela instru¸c˜ao usada foi um ARIMA(3,0,1) com m´edia diferente de zero e AIC= 3966.32, AICc= 3966.43 e BIC= 3993.86. A an´alise dos res´ıduos mostra que a hip´otese de independˆencia dos res´ıduos n˜ao ´e rejeitada (at´e 10 lags), para o teste de Ljung-Box [Ljung e Box (1978)]. Figura 4.11: An´alise de Res´ıduos do Modelo ARIMA. 54 CAP´ ITULO 4. APLICAC¸ ˜ AO Para fornecer uma maior consistˆencia aos resultados, utilizou-se o mesmo teste de Ljung-Box para testar a independˆencia serial na s´erie original, que por sua vez foi rejeitada (note-se que este teste faz uma aproxima¸c˜ao `a s´erie por modelos ARMA). Foi tamb´em utilizado um teste n˜ao param´etrico de Kruskall-Wallis para testar a igualdade de distribui¸c˜oes emp´ıricas entre a vari´avel Dia e vari´avel Velocidade de Vento, que n˜ao foi rejeitada, tendo obtido um valor-p de 0.48. Apesar de os resultados aparentarem aproximar um sinal de ru´ıdo branco (tal como se espera), n˜ao se deve esquecer que os modelos ARIMA usados s˜ao lineares e, existindo padr˜oes n˜ao-lineares, podem ser captados pelo modelo e adulterar as previs˜oes por ele feitas. Por esse motivo, foi explorada uma outra t´ecnica, designada por An´alise Espetral Singular (SSA) [Elsner e Tsonis (2013)], que n˜ao obriga ao conhecimento sobre o modelo param´etrico da s´erie temporal. ´ E composta pelas etapas de decomposi¸c˜ao - por sua vez composta nas etapas de embutimento edecomposi¸c˜ao em valores singulares - e reconstru¸c˜ao - por sua vez composta pelo agrupamento de triplos pr´oprios em´edia da diagonal. Assumindo que a s´erie estudada ´e decomposta na soma de tendˆencia, componentes oscilat´orias e ru´ıdo, este m´etodo s´o obriga `a defini¸c˜ao de dois parˆametros: tamanho da janela L(cujo comprimento ´otimo ´e dado por Lmax =N 2, onde N´e o n´umero de observa¸c˜oes da s´erie; para se conseguir alcan¸car separabilidade suficiente das componente, aconselha-se a escolha de um valor de Lproporcional ao per´ıodo de sazonalidade dos dados [Golyandina et al. (2001)]) a usar e triplos pr´oprios a agrupar. Estes conceitos ser˜ao explicados de seguida. Procedimento SSA: Cada triplo ´e um vetor pr´oprio, um vetor fator e um vetor singular. Na etapa sequencial de SSA, inicialmente faz-se a extra¸c˜ao de tendˆencia usando o primeiro vetor pr´oprio, tal como descrito por Golyandina e Korobeynikov (2014). Em seguida extraem-se as componentes aleat´orias dos res´ıduos. Para tal, devem agrupar-se os triplos pr´oprios com valores singulares pr´oximos, uma vez que a quebra no espetro dos valores pr´oprios permite detetar uma sequˆencia lentamente decrescente de valores singulares produzida por um sinal de ru´ıdo branco. Antes de decidir quais os triplos pr´oprios a agrupar, deve analisar-se a matriz de correla¸c˜oes ponderadas entre as componentes obtidas na separa¸c˜ao de valores pr´oprios (componentes reconstru´ıdas). Deve analisar-se ainda o gr´afico dos vetores pr´oprios sucessivos, agrupando os triplos pr´oprios associados a pol´ıgonos regulares. Note-se que este mecanismo funciona com base no conceito de separabilidade, defendendo que as diferentes componentes da s´erie s˜ao identific´aveis e separ´aveis, permitindo decompor a mesma. Habitualmente toma-se para a extra¸c˜ao de tendˆencia um valor de Lproporcional ao per´ıodo da s´erie, mas como n˜ao se identifica (visualmente) nenhum per´ıodo na s´erie original, considerar-se-´a, nas duas fases do procedimento, que o tamanho da janela ser´a L=728 2= 364 (dias). 4.2. PREVIS ˜ AO - MECANISMOS DE EXPLORAC¸ ˜ AO 55 Figura 4.12: SSA - Valores Singulares. Figura 4.13: SSA - Matriz de Correla¸c˜oes entre as Componentes. Figura 4.14: SSA - Vetores Pr´oprios e Pares de Vetores Pr´oprios. Embora a an´alise gr´afica n˜ao tenha sido trivial, uma vez que apenas o par (4, 5) ´e claramente identific´avel (porque estes dois vetores est˜ao emparelhados na Figura 4.12, apresentam igual comportamento nos vetores pr´oprios da Figura 4.14 e correspondem `as poucas componentes claramente identificadas por uma alta correla¸c˜ao na Figura 4.13), pela an´alise gr´afica destas trˆes componentes de decis˜ao, os triplos escolhidos foram os seguintes: (4, 5) e (7, 8). Observe-se tamb´em que a SSA n˜ao identifica periodicidade na s´erie, uma vez que n˜ao existem pol´ıgonos regulares definidos no gr´afico de Pares de Vetores Pr´oprios da Figura 4.14. Esta informa¸c˜ao ´e ´util para concluir algo que j´a se desconfiava: a s´erie aparenta ser definida por diversas componentes oscilat´orias, cada uma com a sua quota parte de participa¸c˜ao na s´erie original, mas n˜ao parece existir periodicidade em nenhuma das componentes dos triplos selecionados, pelo que a separabilidade da s´erie pode ser feita (usando os triplos pr´oprios e a tendˆencia entretanto extra´ıdos), mas n˜ao se consegue fazer uma interpreta¸c˜ao f´ısica das mesmas. 56 CAP´ ITULO 4. APLICAC¸ ˜ AO Uma vez feito o agrupamento, prosseguiu-se com a reconstru¸c˜ao da s´erie, que se encontra na Figura 4.15, onde se vˆe a s´erie original, a remo¸c˜ao da tendˆencia (a azul) na primeira etapa da SSA sequencial e que centra a s´erie, e a reconstru¸c˜ao da mesma com base nos triplos pr´oprios definidos anteriormente. Figura 4.15: SSA - Reconstru¸c˜ao da S´erie. Com esta informa¸c˜ao, ´e poss´ıvel prever dados com base no algoritmo (de recorrˆencia) descrito no cap´ıtulo 5 da obra de Golyandina et al. (2001). Figura 4.16: SSA - Previs˜ao. 4.3. PREVIS ˜ AO - APLICAC¸ ˜ AO E COMPARAC¸ ˜ AO DE RESULTADOS 57 4.3 Previs˜ao - Aplica¸c˜ao e Compara¸c˜ao de Resultados Para a previs˜ao, tal como referido no Cap´ıtulo 3, foram usadas Redes Neuronais Artificiais, M´aquinas de Suporte Vetorial e ´ Arvores de Regress˜ao (cujas instru¸c˜oes em R s˜ao dadas por nnet, svm e rpart, respetivamente). Relembrando, a valida¸c˜ao cruzada ´e usada para avaliar a capacidade de predi¸c˜ao/previs˜ao do modelo com uns determinados parˆametros (o k-fold ´e o m´etodo de parti¸c˜ao dos dados e ´e independente do m´etodo que se usa para a separa¸c˜ao inicial; ´e usado para que haja reamostragem e o treino n˜ao seja feito sempre sobre os mesmos dados). Usaram-se 10 folds para a valida¸c˜ao, o que significa que o algoritmo na valida¸c˜ao cruzada ´e treinado com 90% do treino e testado nos restantes 10% do treino. O erro de previs˜ao final do modelo, nesta etapa, ´e dado pelo MAPE no conjunto de treino (embora tamb´em tenha sido calculado o erro MAD). Com a informa¸c˜ao que se obt´em do passo anterior, determinam-se os parˆametros que obtˆem um menor erro no treino (isto para cada um dos m´etodos estudados). Nesta etapa, j´a se disp˜oe dos valores de desempenho de cada modelo. Comparando-os, existe um modelo vencedor e ´e esse modelo que vai ser usado para fazer previs˜oes no conjunto de teste. Os conjuntos de treino e teste iniciais foram definidos de duas formas: •Aleat´oria: separa¸c˜ao aleat´oria da amostra em dois conjuntos: treino - 70% e teste - 30%. •Sequencial: separa¸c˜ao da amostra em dois conjuntos: treino - primeiros 70% das observa¸c˜oes da base de dados - e teste - restantes 30% das observa¸c˜oes da base de dados. Uma vez que esta separa¸c˜ao ´e organizada no tempo, ´e poss´ıvel comparar os modelos acima vistos com o modelo ARIMA, usando a t´ecnica de Sliding Window. Far-se-´a de seguida a an´alise dos resultados obtidos. An´alise dos resultados da separa¸c˜ao inicial aleat´oria: Os resultados da valida¸c˜ao cruzada para cada modelo foram os seguintes: •´ Arvores de regress˜ao: Figura 4.17: Previs˜oes Num´ericas por ´ Arvores de Regress˜ao. 64 CAP´ ITULO 5. MODELO GENERALIZADO DE SPARRE ANDERSEN Assume-se que o processo do n´umero de indemniza¸c˜oes, dado por {Nt:t= 0,1, ...}, ´e um processo de renovamento modificado em tempo discreto, com tempos entre indemniza¸c˜oes positivos e independentes, onde W1´e a dura¸c˜ao do tempo 0 at´e ao instante da primeira indemniza¸c˜ao e Wio tempo entre as (i-´esima-1) e i-´esima indemniza¸c˜ao. H´a que notar que (W1,W2,...) forma uma sequˆencia de vari´aveis aleat´orias independentes e identicamente distribu´ıdas com fun¸c˜ao de probabilidade dada por aj=P(Wi=j)j= 1,2, ..., na(na<∞) e correspondente fun¸c˜ao de sobrevivˆencia Aj=P(Wi> j) = 1 − j X k=1 ak. Assume-se tamb´em que a distribui¸c˜ao do tempo entre indemniza¸c˜oes tem suporte finito. No contexto da an´alise de risco feita numa seguradora, segundo o modelo cl´assico, os tempos entre indemniza¸c˜oes formam uma sequˆencia de vari´aveis aleat´orias independentes com distribui¸c˜ao exponencial de parˆametro λ. Assim, o n´umero de indemniza¸c˜oes segue um processo de Poisson, o que, devido `a perda de mem´oria que caracteriza a a exponencial, faz com que as vari´aveis “tempo at´e `a primeira indemniza¸c˜ao” (W1) e “instante da primeira indemniza¸c˜ao” (T1) tenham a mesma distribui¸c˜ao. No caso do GMSA, a forma como este processo ´e conduzido ´e geral, na medida em que a distribui¸c˜ao dos tempos entre indemniza¸c˜oes ´e arbitr´aria (e portanto a distribui¸c˜ao de W1n˜ao ser´a necessariamente a mesma de T1, a n˜ao ser que tenha ocorrido alguma indemniza¸c˜ao no instante 0). Ora, como W1ser´a tratado como um caso `a parte, uma vez que n˜ao existe informa¸c˜ao sobre o que se passou antes da ocorrˆencia da primeira indemniza¸c˜ao, h´a duas considera¸c˜oes feitas habitualmente: •assume-se que ocorreu uma indemniza¸c˜ao antes de 0, pelo que W1,W2,... passam a ter a mesma distribui¸c˜ao •assume-se o modelo mais geral, em que W1segue uma qualquer fun¸c˜ao de probabilidade rj=P(W1=j)j= 1,2, ..., nr(nr<∞) e correspondente fun¸c˜ao de sobrevivˆencia Rj=P(W1> j) = 1 − j P k=1 rk Assim, o principal processo de risco analisado ´e o seguinte: Ut=u+ t−1 X i=0 pi− Nt X i=1 Yi onde Utrepresenta o montante do super´avit da seguradora no instante t, ou seja, a reserva de risco de uma carteira no instante t. Este processo ´e, na verdade, um balan¸co de contas entre o que entra de pr´emios na seguradora e o que sai para pagamento de indemniza¸c˜oes, 5.1. QUANTIDADES DE RU´ INA E LIMIAR DE PR´ EMIO 65 n˜ao esquecendo a reserva inicial dada por u=U0∈N. Da forma como se define, Utrepresenta o montante de super´avit no final do intervalo (t−1, t]. Em rela¸c˜ao a este intervalo em particular, assume-se que os pr´emios s˜ao recebidos em (t−1)+ e as indemniza¸c˜oes pagas em t−, garantindo um balan¸co real em t. ´ E agora adicionada ao super´avit da seguradora uma quantia Υ ∈Z+que afeta o pr´emio recebido num certo instante e que funciona como um limiar entre um pr´emio constante e um pr´emio aleat´orio (isto porque, a partir de um determinado montante de reserva, a companhia j´a tem seguran¸ca financeira suficiente para flexibilizar os pr´emios dos seus clientes). Assim, define-se que um pr´emio aleat´orio pt=cse Ut<Υ Xtse Ut≥Υ com c∈Z+como o pr´emio sem encargos e Xtcomo um pr´emio aleat´orio recebido em t, com di=P(Xt=i), i =c1, c1+ 1, ..., c2, c2 X i=c1 di= 1. c1, c2s˜ao, respetivamente, os valores m´ınimo e m´aximo do suporte da distribui¸c˜ao de Xt, com c1, c2∈ {0,1, ..., c}ec1≤c2. Assim, a quantia c−ptpode ser interpretada como a participa¸c˜ao de resultados, muitas vezes referida no mundo segurador. Sejam {Y1, Y2, Y3...}os montantes individuais das indemniza¸c˜oes que a seguradora tem de pagar - vari´aveis aleat´orias positivas i.i.d. com fun¸c˜ao densidade de probabilidade αj=P(Yi=j)j= 1,2, ..., mα e correspondente fun¸c˜ao de sobrevivˆencia Λj= 1 − j X k=1 αk(mα≤ ∞). Consideram-se as seguintes quantidades de ru´ına: •oinstante de ru´ına, escrito como T=min{t∈Z+|Ut<0} eT=∞se Ut≥0∀t∈Z+ •se a ru´ına ocorrer, o d´efice na ru´ına ´e dado por |UT| •se a ru´ına ocorrer, o super´avit imediatamente anterior `a ru´ına ´e dado por UT−=UT−1+pT−1 66 CAP´ ITULO 5. MODELO GENERALIZADO DE SPARRE ANDERSEN Portanto T=∞se mα≤min{c, Υ + c1}. Por outro lado, se mα> min{c, Υ + c1}, ent˜ao |UT|∈ {1,2, ..., mα−min{c, Υ + c1}} e UT−∈ {min{c, Υ + c1}, min{c, Υ + c1}+ 1, ..., mα−1}}. O que a nota¸c˜ao acima definida mostra ´e que, se o valor das indemniza¸c˜oes a pagar aos segurados ´e superior `a entrada de dinheiro que existe na seguradora, proveniente do pagamento de pr´emios, a ru´ına ´e certa e o d´efice da ru´ına pode ir desde uma unidade monet´aria at´e `a diferen¸ca entre o valor que a companhia teve de pagar em indemniza¸c˜oes e o montante de que dispunha at´e esse acontecimento. Do mesmo modo, o super´avit imediatamente anterior `a ru´ına pode ir desde o valor m´aximo de que a seguradora dispunha antes da indemiza¸c˜ao (caso mα=min{c, Υ + c1}) at´e a uma unidade monet´aria a menos que a necess´aria para pagar essa mesma indemniza¸c˜ao. Uma vez definidas as quantidades de ru´ına a modelar, devem introduzir-se as probabilidades conjuntas que lhes est˜ao associadas e onde se pretende chegar para que a modela¸c˜ao seja poss´ıvel. Sejam: •ωn,i(u) = P(T=n, UT−=i, |U0=u) •φn,j(u) = P(T=n, |UT|=j|U0=u) •ψn,i,j(u) = P(T=n, UT−=i, |UT|=j|U0=u) Mais `a frente, ser´a visto o modo de estimar as probabilidades associadas com o uso das vari´aveis aleat´orias |UT|eUT−. 5.2 Formula¸c˜ao do Modelo Daqui em diante, seja W=Wiarbitr´ario com i= 2,3, ... e τj=P(W > j |W > j −1) = Aj Aj−1 (τ1=A1eτna−1= 0) a probabilidade do tempo entre indemniza¸c˜oes ser superior a junidades de tempo, sabendo que at´e ao instante de an´alise anterior, a indemniza¸c˜ao n˜ao tinha ocorrido. Considera-se ainda S=         0τ10··· 0 0 0 τ2··· 0 . . .. . ........ . . 0 0 ...0τna−1 0 0 ...0 0         ,s=       1−τ1 1−τ2 . . . 1−τna−1 1        ee1=1,0, ..., 0 tais que aj=P(Wi=j) = e1Sj−1s. Esta defini¸c˜ao foi dada e demonstrada por Wu e Li (2008). 5.2. FORMULAC¸ ˜ AO DO MODELO 67 O objetivo interm´edio deste modelo ´e a constru¸c˜ao do processo bivariado {(Ut, Lt) : t=k, k + 1, ...}. em que Utrepresenta - tal como visto anteriormente - o super´avit da seguradora no instante teLtdenota um contador de tempo em tque mede o “tempo que falta” at´e `a pr´oxima indemniza¸c˜ao. Note-se ainda que a componente Ucorresponde ao n´ıvel do processo, ao passo que a componente Lcorresponde `a fase do processo, que segue uma rela¸c˜ao Markoviana: (Ut+1, Lt+1) = (Ut+pt, Lt+ 1) se n˜ao existe indemniza¸c˜ao em (t+ 1)+ (Ut+pt−Y, 1) se existe uma indemniza¸c˜ao de Yem (t+ 1)+ ´ E agora poss´ıvel analisar a matriz que cont´em as probabilidades de transi¸c˜ao associadas a esta cadeia de Markov (com um espa¸co de estados dado por ∆ = Zx{1,2, ..., na}): ··· −1 0 1 ··· Υ−2 Υ −1 Υ ···                                   . . ..... . .. . .. . .··· . . .. . .. . .. . . −1··· BcBc−1Bc−2··· Bc−Υ+1 Bc−ΥBc−Υ−1··· 0··· Bc+1 BcBc−1··· Bc−Υ+2 Bc−Υ+1 Bc−Υ··· 1··· Bc+2 Bc+1 Bc··· Bc−Υ+3 Bc−Υ+2 Bc−Υ+1 ··· . . .··· . . .. . .. . ..... . .. . .. . .. . . Υ−2··· BΥ+c−1BΥ+c−2BΥ+c−3··· BcBc−1Bc−2··· Υ−1··· AΥ+1 AΥAΥ−1··· A2A1A0··· Υ··· AΥ+2 AΥ+1 AΥ··· A3A2A1··· Υ+1 ··· AΥ+3 AΥ+2 AΥ+1 ··· A4A3A2··· . . .··· . . .. . .. . .··· . . .. . .. . .... onde cada elemento da matriz ´e uma de duas matrizes por blocos, definidas como Bi=   Onaif i∈Z− Sse i= 0 (se1)αise i∈Z+ eAi= c2 X j=c1 djBi+j Note-se que este processo ´e uma cadeia de Markov dupla, cuja matriz Papresenta blocos finitos de dimens˜ao naxna. Como cada elemento da matriz corresponde a outra matriz, a Figura 5.1 seguinte tenta ilustrar a estrutura da matriz, para que haja uma melhor interpreta¸c˜ao da mesma. 68 CAP´ ITULO 5. MODELO GENERALIZADO DE SPARRE ANDERSEN Figura 5.1: Estrutura Matricial de Transi¸c˜oes no GMSA. Considere-se agora um espa¸co de estados onde ∆1=Nx{1,2, ..., na}e ∆2=Z−x{1,2, ..., na} s˜ao tais que ∆1∩∆2=∅. Tomem-se agora duas matrizes de transi¸c˜ao (sub-matrizes de P) C: ∆1→∆1eD: ∆1→∆2 que mapeiam, respetivamente, estados de n˜ao ru´ına em ∆1e estados de ru´ına em ∆2(C´e o quadrante inferior direito de PeD´e o quadrante inferior direito de Phorizontalmente invertida). Para o c´alculo das quantidades relacionadas com a ru´ına e vistas anteriormente, ter˜ao de ser acrescentadas algumas defini¸c˜oes, necess´arias mais `a frente. Sejam: •zt=min{i∈ {1,2, ..., t} | u+c(i−1)}se u < T 0 se u≥T •dm,n =   c2 P j=c1 dj,1dm−j,n−1se m=nc1, nc1+ 1, ..., nc2 0 caso contr´ario lembrando que dm,0=δm,0(onde δ´e a fun¸c˜ao de Kronecker) dm,1=dmse m=nc1, nc1+ 1, ..., nc2 0 caso contr´ario Como W1=k,b(k)´e o vector-linha com as probabilidades iniciais de estar nos estados de ∆1e ´e dada por (k−zk)c2 X m=(k−zk)c1 dm,k−zk(αu+czk+me1, αu+czk+m−1e1, ..., α2e1, α1e1,0,0, ...). 5.2. FORMULAC¸ ˜ AO DO MODELO 69 Agora s˜ao definidos dois vectores-coluna fundamentais: •g(k) n= (g(k) n,0,g(k) n,1,g(k) n,2, ...) = b(k)Cnn∈N, correspondente `a probabilidade de transi¸c˜ao de um estado de n˜ao ru´ına para outro de n˜ao ru´ına em kunidades de tempo; •h(k) n= (h(k) n,−1,h(k) n,−2,h(k) n,−3, ...) = b(k)Cn−1D n ∈Z+, que corresponde `a probabilidade de transi¸c˜ao de um estado de n˜ao ru´ına para um estado de ru´ına em kunidades de tempo. Lembrando que h(k) n,−j= (φ(k) n,j(u),0, ..., 0) com n∈Z+ φ(k) n,j(u) = P(T=k+n, |UT|=j|Uk∈Ωk) j= 1,2, ..., mα−min{c, Υ + c1}, Ωk={0,1, ..., u +czk+c2(k−zk)} obt´em-se φ(k) n,j(u) = h(k) n,−jeT 1. Aplicando o mesmo racioc´ınio, obt´em-se uma representa¸c˜ao para ψ(k) n,i,j(u) = P(T=k+n, UT=i, |UT|=j|Uk∈Ωk) de onde segue imediatamente que ψ(k) n,i,j(u) = (g(k) n−1,i−pk+n−1s)αi+j. Uma vez que ´e utilizado um limiar Υ e que as quantidades de ru´ına podem tomar valores muito particulares (como se ver´a), a ´ultima probabilidade pode ser calculada com base nas express˜oes j´a vistas, dada pelo seguinte corol´ario [Drekic e Mera (2011)]: Corol´ario para calcular g(k) n−1,i−pk+n−1 •se min{c, Υ + c1}=c, ent˜ao              g(k) n−1,i−cse i=c, c + 1, ..., Υ + c1−1 g(k) n−1,i−c+ min{c2,i−Υ} P j=c1 djg(k) n−1,i−jse i= Υ + c1, ..., min{mα,Υ + c}−1 c2 P j=c1 djg(k) n−1,i−jse i= Υ + c, Υ + c+ 1, ..., mα−1 70 CAP´ ITULO 5. MODELO GENERALIZADO DE SPARRE ANDERSEN •se min{c, Υ + c1}= Υ + c1, ent˜ao                  min{c2,i−Υ} P j=c1 djg(k) n−1,i−jse i= Υ + c1,Υ + c1+ 1, ..., c −1 g(k) n−1,i−c+ min{c2,i−Υ} P j=c1 djg(k) n−1,i−jse i=c, c + 1, ..., min{mα,Υ + c}−1 c2 P j=c1 djg(k) n−1,i−cse i= Υ + c, Υ + c+ 1, ..., mα−1 Lembrando que W1=k, k ∈ {1,2, ..., nr}, a Lei total das Probabilidades leva a que se conclua que •φn,j(u) = nr P k=1 rkφ(k) n−k,j(u) •ψn,i,j(u) = nr P k=1 rkψ(k) n−k,i,j(u) n=nr+ 1, nr+ 2, ... No entanto, pode acontecer que T=npara n= 1,2, ..., nr, caso que n˜ao se deve esquecer. H´a duas explica¸c˜oes poss´ıveis para que aconte¸ca: 1. a primeira indemniza¸c˜ao ocorre em n−para um limite de super´avit de u+czn+m (onde n−1 P i=zn Xi=m); 2. a primeira indemniza¸c˜ao ocorre num instante k−k∈ {1,2, ..., n−1}e a ru´ına ocorre n−kunidades de tempo depois. Combinando a informa¸c˜ao anterior com a que resulta para os casos em que n=nr+ 1, nr+ 2, ... (n∈Z+), definem-se agora as express˜oes gerais para as quantidades de ru´ına que se pretendem calcular (tal como descrito no in´ıcio da explica¸c˜ao do modelo): •φn,j(u) = min{n−1,nr} P k=1 rkφ(k) n−k,j(u) + rn (n−zn)c2 P m=(n−zn)c1 dm,n−znαu+czn+m+j; •ψn,i,j(u) = min{n−1,nr} P k=1 rkψ(k) n−k,i,j(u) + rndi−u−czn,n−znαi+j. Notando que Λu+cn+mα−i= 0 se mα=∞, •a fun¸c˜ao de massa bivariada ωn,i(u) ´e escrita como min{n−1,nr} X k=1 rk mα−i X j=1 ψ(k) n−k,i,j(u) + rndi−u−czn,n−zn mα−i X j=1 αi+j! = Λi  min{n−1,nr} X k=1 rk(g(k) n−k−1,i−pn−1s) + rndi−u−czn,n−zn . 5.2. FORMULAC¸ ˜ AO DO MODELO 71 Concluindo a explica¸c˜ao, ´e de referir que, para que seja poss´ıvel analisar todo o processo, e uma vez que as express˜oes imediatamente acima resultam da soma das j´a obtidas para todos os valores de kunidades de tempo a considerar (e que depender˜ao obviamente de cada situa¸c˜ao que se pretenda analisar), ser´a necess´ario somar todos os valores que se pretendam estudar para o d´efice na ru´ına esuper´avit imediatamente anterior `a ru´ına, obtendo assim a fun¸c˜ao de distribui¸c˜ao trivariada Ψn,x,y(u), definida como P(T≤n, UT−≤x, |UT|≤ y|U0=u) = n X l=1 x X i=min{c,Z+c1} y X j=1 ψl,i,j(u). Este modelo foi inicialmente pensado como principal foco de trabalho, tendo sido iniciada a sua implementa¸c˜ao. No entanto, devido `as dificuldades enfrentadas nos processos de modela¸c˜ao e previs˜ao de velocidade m´axima de ventos, optou-se pela sua defini¸c˜ao te´orica, que, apesar de n˜ao ser original, salienta a capacidade de generaliza¸c˜ao associada `a modela¸c˜ao de um processo t˜ao complexo quanto o de montantes geridos por uma Seguradora, que se constitui como um dos processos de gest˜ao mais lucrativos que se conhece. Ilustra¸c˜ao: Apresentam-se de seguida alguns resultados obtidos at´e ao momento pela implementa¸c˜ao parcial deste modelo em R (com base em exemplos definidos em Drekic e Mera (2011)). Para o caso apresentado, assumem-se os seguintes parˆametros: •f.d.p. geom´etrica truncada aj=(2 9)9 11)j−1se j= 1,2, ..., na−1 9 11)na−1se j=na para na= 10 •f.d.p. pareto αj=1 + j−1 30 )−4−1 + j 30)−4para j∈Z+ema= 100 •di=5 i2 5)i3 5)5−i, para i= 0, ..., c •c= 5, u= 2, Υ = 50 72 CAP´ ITULO 5. MODELO GENERALIZADO DE SPARRE ANDERSEN Apresenta-se de seguida uma parte do output da matriz P, definida por blocos: Figura 5.2: P definida por Blocos e Extra¸c˜ao de C e D - Exemplo Implementado. FCUP 73 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Cap´ıtulo 6 Trabalho Futuro Janssen e Manca (1997), que s˜ao os principais impulsionadores da aplica¸c˜ao de modelos semi-markovianos `as mais variadas ´areas, e que detˆem diversas publica¸c˜oes sobre o tema, escrevem que o “Sob o ponto de vista computacional” (...) “´e claro que este modelo n˜ao pode ser utilizado como um modelo de simula¸c˜ao sem ser na presen¸ca de uma boa m´aquina”1. Efetivamente, estes s˜ao modelos conceptualmente e computacionalmente “pesados”, respetivamente devido `a quantidade de informa¸c˜ao de que disp˜oem e devido ao facto de terem, por norma, processos definidos por recorrˆencia (como se pode ver na equa¸c˜ao de evolu¸c˜ao do processo semi-markoviano de primeira ordem). Os modelos semi-markovianos s˜ao uma generaliza¸c˜ao dos modelos de markov, motivo pelo qual d˜ao uma maior liberdade na abordagem em qualquer tipo de problema que comporte modela¸c˜ao de processos que dependam do tempo, filas de espera, etc. Como trabalho futuro seria interessante executar o mesmo tipo de simula¸c˜ao feita na modela¸c˜ao tendo em conta a recorrˆencia temporal no processo de segunda ordem e avaliando as suas diferen¸cas comportamentais em compara¸c˜ao com as restantes simula¸c˜oes. Seria tamb´em poss´ıvel a considera¸c˜ao de vari´aveis explicativas no processo de modela¸c˜ao das velocidades de vento (por exemplo, as que se declararam significativas nos modelos de regress˜ao), usando-as na estima¸c˜ao de P e F. Relativamente `a modela¸c˜ao da fun¸c˜ao densidade de probabilidade, e se se dispusesse de uma amostra temporal mais alargada, poder-se-ia recorrer a outros m´etodos de estima¸c˜ao da mesma, param´etricos ou n˜ao param´etricos, ou at´e uma conjuga¸c˜ao dos dois (ver, por exemplo, McNeil (1999)). Do ponto de vista da previs˜ao, uma an´alise mais aprofundada dos modelos ARIMA poderia sugerir melhores resultados a partir de um estudo temporal diferente dos mesmos dados (semanal, quinzenal, trimestral, mensal, avalia¸c˜ao por esta¸c˜oes do ano, etc.). Note-se ainda que, para o tipo de problema considerado, dois anos de dados podem n˜ao ser suficientemente informativos sobre o processo, compromentendo a efic´acia da previs˜ao. Para al´em do trabalho efetuado, poder-se-iam considerar modelos n˜ao lineares para a previs˜ao. 1tradu¸c˜ao livre da autora. 80 BIBLIOGRAFIA Ross SM (2014). Introduction to probability models. Academic Press. (p´ag. 19) Sansom J, Thomson P, et al. (2001). “Fitting hidden semi-Markov models to breakpoint rainfall data.” Journal of Applied Probability,38, 142–157. (p´ag. 19) Suykens JA, Vandewalle JP, de Moor BL (2012). Artificial neural networks for modelling and control of non-linear systems. Springer Science & Business Media. (p´ag. 33) Torfs P, Brauer C (2014). “A (very) short introduction to R.” Hydrology and Quantitative Water Management Group Wageningen University. (p´ag. 87) Torgo L (2010). Data mining with R: learning with case studies. Chapman & Hall/CRC. (p´ag. 38), (p´ag. 87) Vapnik VN (1995). The nature of statistical learning theory. Springer-Verlag New York, Inc. (p´ag. 34) Watson GS (1964). “Smooth regression analysis.” Sankhy¯a: The Indian Journal of Statistics, Series A,26, 359–372. (p´ag. 34) Witten IH, Frank E (2005). Data Mining: Practical machine learning tools and techniques. Morgan Kaufmann. (p´ag. 35) Wu X, Li S (2008). On a discrete-time Sparre Andersen model with phase-type claims. Centre for Actuarial Studies: the research paper series, No. 169, University of Melbourne. (p´ag. 66) FCUP 81 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Anexos FCUP 83 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Anexo A Demonstra¸c˜ao. (referente `a pagina 21). Seja Pn(t) = P(N(t) = n). Considere-se P0(t+h) = P(N(t+h)=0)=P(N(t) = 0, N(t+h)−N(t) = 0) =P(N(t) = 0)P(N(t+h)−N(t) = 0) = P0(t)(1 −λh +o(h)). Logo, obt´em-se P0 0(t) = lim h→0 P0(t+h)−P0(t) h= lim h→0−λP0(t) + o(h) h=−λP0(t). A solu¸c˜ao desta equa¸c˜ao ´e P0(t) = Ce−λt. Usando a condi¸c˜ao inicial P0(0) = 1, obt´em-se C= 1 e conclui-se que P0(t) = e−λt. Para n > 0, obt´em-se: Pn(t+h) = Pn(t)P(N(t+h)−N(t) = 0) + Pn−1(t)P(N(t+h)−N(t) = 1)+ + n X k=2 Pn−k(t)P(N(t+h)−N(t) = k) =Pn(t)(1 −λh +o(h)) + Pn−1(t)(λh +o(h)) + n X k=2 Pn−k(t)o(h) Logo P0 n(t) = lim h→0 Pn(t+h)−Pn(t) h=−λPn(t) + λPn−1(t)+ + lim h→0Pn(t)o(h) h+Pn−1(t)o(h) h+ n X k=2 Pn−k(t)o(h) h=−λPn(t) + λPn−1(t). Reescrevendo a equa¸c˜ao, tem-se eλt(P0 n(t) + λPn(t)) = λeλtPn−1(t) ou d dt(eλtPn(t)) = λeλtPn−1(t). Assim, d dt(eλtPn(t)) = λntn−1 (n−1)! ou eλtPn(t) = (λt)n n!+C. Pela condi¸c˜ao inicial C= 0, logo Pn(t) = e−λt (λt)n n!. 84 Informa¸c˜ao de apoio `a decis˜ao Aqui est˜ao inclu´ıdas algumas an´alises gr´aficas ou num´ericas de outputs das diferentes instru¸c˜oes usadas em R ao longo desta tese. Esta informa¸c˜ao ´e referida de modo contextualizado, referente ao Cap´ıtulo de Aplica¸c˜ao: •CAP´ ITULO 4 - APLICAC¸ ˜ AO: –Resultado da aplica¸c˜ao do m´etodo dos m´ınimos generalizados no caso hor´ario da primeira etapa de an´alise da base de dados: Generalized least squares fit by REML Model: VMax ~ I(A) + I(M) + I(H) + I(E) + TMin + HMax + RMin + DV + PT Data: BaseDados AIC BIC logLik 92031.63 92132.81 -46002.81 Coefficients: Value Std.Error t-value p-value (Intercept) -193.63100 94.41517 -2.05085 0.0403 I(A) 0.10226 0.04686 2.18209 0.0291 I(M) 0.02959 0.01328 2.22820 0.0259 I(H) 0.00537 0.00356 1.51173 0.1306 I(E)Outono 0.49885 0.12423 4.01566 0.0001 I(E)Primavera -0.79088 0.08837 -8.94990 0.0000 I(E)Ver~ao -0.30683 0.12943 -2.37060 0.0178 TMin -0.18092 0.00716 -25.28111 0.0000 HMax -0.03994 0.00152 -26.30061 0.0000 RMin 0.00220 0.00015 14.27522 0.0000 DV -0.01349 0.00022 -61.81387 0.0000 PT 0.46016 0.02432 18.92248 0.0000 –Resultado da aplica¸c˜ao da regress˜ao linear no caso di´ario da primeira etapa de an´alise da base de dados: Call: lm(formula = VMax ~ I(E) + TMed + HMed + DV + PT,data = d) Residuals: Min 1Q Median 3Q Max -9.3278 -2.3320 -0.3839 2.0254 14.5806 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 18.825446 0.679681 27.697 < 2e-16 *** I(d$Estacao)Outono 1.269779 0.426506 2.977 0.003007 ** I(d$Estacao)Primavera 0.237417 0.460728 0.515 0.606496 I(d$Estacao)Ver~ao 1.818327 0.612668 2.968 0.003098 ** d$Tmed -0.403030 0.034928 -11.539 < 2e-16 *** d$HMed -0.031680 0.008215 -3.856 0.000125 *** d$DV -0.014155 0.001664 -8.505 < 2e-16 *** d$PT 0.169784 0.071727 2.367 0.018192 * 85 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.488 on 720 degrees of freedom Multiple R-squared: 0.4042,Adjusted R-squared: 0.3984 F-statistic: 69.77 on 7 and 720 DF, p-value: < 2.2e-16 –Resultado da aplica¸c˜ao das fun¸c˜oes de modela¸c˜ao de variˆancia - melhor modelo obtido: Generalized least squares fit by REML Model: VMax ~ I(A) + I(M) + I(Estacao) + TMin + HMax + RMin + DV + PT Data: BaseDados AIC BIC logLik 88984.99 89109.52 -44476.5 Variance function: Structure: Exponential of variance covariate, different strata Formula: ~TMin | E Parameter estimates: Ver~ao Outono Inverno Primavera -0.04943755 -0.03987934 -0.05740226 -0.05254327 Coefficients: Value Std.Error t-value p-value (Intercept) -620.2921 77.89434 -7.96325 0.0000 I(A) 0.3139 0.03866 8.11979 0.0000 I(M) 0.0840 0.01374 6.11193 0.0000 I(Estacao)Outono 0.5408 0.12735 4.24620 0.0000 I(Estacao)Primavera -0.7176 0.08818 -8.13741 0.0000 I(Estacao)Ver~ao -0.3545 0.11579 -3.06165 0.0022 TMin -0.1638 0.00631 -25.94540 0.0000 HMax -0.0491 0.00139 -35.30828 0.0000 RMin 0.0018 0.00010 17.51989 0.0000 DV -0.0106 0.00019 -54.70199 0.0000 PT 0.4684 0.02762 16.95941 0.0000 Correlation: (Intr) I(A) I(M) I(Es)O I(Es)P I(Es)V TMin HMax RMin DV I(A) -1.000 I(M) -0.251 0.250 I(Estacao)O -0.057 0.057 -0.678 I(Estacao)P 0.150 -0.150 -0.236 0.608 I(Estacao)V -0.132 0.133 -0.395 0.763 0.779 TMin -0.036 0.035 -0.091 -0.245 -0.536 -0.631 HMax 0.049 -0.050 0.025 -0.205 -0.221 -0.228 0.215 RMin -0.037 0.037 0.079 0.058 0.098 0.164 -0.437 0.258 DV -0.027 0.028 -0.017 0.062 0.029 0.015 -0.042 -0.280 0.015 PT -0.008 0.008 -0.046 0.033 0.032 0.027 0.054 -0.120 -0.040 0.037 Standardized residuals: Min Q1 Med Q3 Max -3.9028458 -0.6491694 -0.1364726 0.5777264 6.6580570 Residual standard error: 6.743717 86 Seguidamente s˜ao apresentados os resultados da previs˜ao temporal (hor´aria) por uso do modelo linear ARMA (com base na transforma¸c˜ao de dados e determina¸c˜ao da componente sazonal do modelo) e por uso de redes neuronais artificiais, com 1 camada escondida (e divis˜ao de dados em 70% para teste, 15% para treino e 15% para valida¸c˜ao). Devido `a alta correla¸c˜ao residual, optou-se pela coloca¸c˜ao destes resultados apenas para visualiza¸c˜ao. Note-se que a componente sazonal parece ter um bom ajuste aos dados, mas foi estimada com base numa s´erie harm´onica (combina¸c˜ao linear de senos e cossenos) por regress˜ao linear. Ora, j´a foi visto que a regress˜ao linear n˜ao deve ser considerada nestes dados, uma vez que n˜ao apresentam independˆencia. Figura 1: Modelo Arima em Dados Hor´arios. Figura 2: Redes Neuronais Artificiais em Dados Hor´arios. FCUP 87 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Anexo B Ferramentas de trabalho Para a explora¸c˜ao da base de dados e implementa¸c˜ao dos modelos j´a referidos, foram utilizadas duas ferramentas de programa¸c˜ao: R e MATLAB. OR[R Core Team and others (2012)] ´e um ambiente de desenvolvimento integrado orientado a objetos, direcionado sobretudo para an´alise e manipula¸c˜ao de dados. Entre as suas principais vantagens, destacam-se o facto de ser uma aplica¸c˜ao de distribui¸c˜ao gratuita e de c´odigo p´ublico, existindo vers˜oes j´a compiladas para os principais sistemas operativos, que, na sua maioria, s˜ao fornecidas pela comunidade de utilizadores. ´ E particularmente ´util para lidar com grandes conjuntos de dados [Torgo (2010)] e apresenta bons tempos de execu¸c˜ao nessas circunstˆancias [Matloff (2011)]. No entanto, existe uma desvantagem de relevo: a existˆencia de packages pr´e-implementados ou definidos pelos utilizadores - e dispon´ıveis para os restantes - obriga ao conhecimento da teoria que est´a por tr´as das implementa¸c˜oes, para que se possa fazer pleno uso desses packages. O que se pretende dizer com isto ´e que o uso de pr´e-implementa¸c˜oes n˜ao deve ser feito de forma despreocupada e sem tentar perceber qual o racioc´ınio que lhes est´a associado, quando este n˜ao ´e ´obvio. O c´odigo fonte para o ambiente de software ´e escrito principalmente em C, FORTRAN e R. Para uma programa¸c˜ao mais intuitiva neste ambiente, optou-se pela utiliza¸c˜ao do RStudio [Torfs e Brauer (2014)], que ´e um ambiente de desenvolvimento integrado (IDE) para o R e que, por esse motivo, apresenta caracter´ısticas e ferramentas de apoio ao desenvolvimento de software com o objetivo de agilizar este processo (como autocomplete de comandos, indica¸c˜oes gr´aficas de fecho e abertura de parˆentesis, indica¸c˜oes de nomes internos de fun¸c˜oes, disponibiliza¸c˜ao f´acil do c´odigo das intru¸c˜oes, etc.). O MATLAB (diminutivo de MATrix LABoratory), ´e um software interativo particularmente indicado para o c´alculo num´erico e usado desde a d´ecada de 70 [Moler et al. (1980), Haigh (2008)], que integra an´alise num´erica, c´alculo matricial, m´etodos de processamento de sinais e constru¸c˜ao gr´afica. Esta ferramenta foi utilizada como valida¸c˜ao de alguns resultados obtidos, tendo sido usadas as toolbox de Redes Neuronais e a aplica¸c˜ao Semi-Markov. Esta ´ultima foi utilizada para poder comparar os resultados provenientes do algoritmo implementado pela autora desta tese e os resultados provenientes dessa toolbox, vistos no gr´afico da esquerda da Figura 4.5. 88 FCUP 89 Modela¸c˜ao e Previs˜ao de Velocidade de Ventos Anexo C Categoriza¸c˜ao da velocidade de vento Existe mais do que uma escala de categoriza¸c˜ao de vento, entre as quais a escala Beaufort, que ´e um sistema que relaciona a velocidade do vento com as condi¸c˜oes observadas no mar ou em terra, isto ´e, com os efeitos da mesma sobre estes dois meios f´ısicos. Os crit´erios est˜ao descritos nas Tabelas 1e2e foram baseados na informa¸c˜ao que se pode encontrar no site da Organiza¸c˜ao Mundial de Meteorologia https://www.wmo.int/pages/ index_en.html 1. Note-se que na base de dados em estudo, na esta¸c˜ao meteorol´ogica de Settala, cuja localiza¸c˜ao geogr´afica se mostra em seguida, o valor m´aximo de velocidade de vento ´e de 25 m/s, n˜ao se justificando a considera¸c˜ao de intensidades de vento elevadas. No entanto, ponderou-se a discretiza¸c˜ao (separa¸c˜ao de valores por estado) aliada a um crit´erio climatol´ogico, tamb´em descrito abaixo, na Tabela 3. Intensidade Descri¸c˜ao mph m/s 0 Calmo <1<0.3 1 Aragem 1 - 3 0.3 - 1.5 2 Brisa Leve 4 - 7 1.6 - 3.3 3 Brisa Fraca 8 - 12 3.4 - 5.4 4 Brisa Moderada 13 - 18 5.5 - 7.9 5 Brisa Forte 19 - 24 8 - 10.7 6 Vento Suave 25 - 31 10.8 - 13.8 7 Vento Forte 32 - 38 13.9 - 17.1 8 Ventania 39 - 46 17.2 - 20.7 9 Ventania Forte 47 - 54 20.8 - 24 10 Tempestade 55 - 63 24.5 - 28.4 11 Tempestade Violenta 64 - 74 28.5 - 32.6 12 Furac˜ao ≥74 ≥32.7 Tabela 1: Escala de Beaufort - Intensidades de vento. 1´ultima consulta realizada em 30 de Agosto de 2015. 96 HMax: Humidade Relativa M´axima Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 8.7 •M´edia amostral: 66.93 •M´aximo amostral: 100.0 •1ºQuartil: 54.3 •2ºQuartil: 69.3 •3ºQuartil: 80.7 RMed: Radia¸c˜ao Global M´edia Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 0.0 •M´edia amostral: 149.1 •M´aximo amostral: 988.0 •1ºQuartil: 0.0 •2ºQuartil: 4.4 •3ºQuartil: 220 RMin: Radia¸c˜ao Global M´ınima Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 0.0 •M´edia amostral: 103.5 •M´aximo amostral: 958.9 •1ºQuartil: 0.0 •2ºQuartil: 0.0 •3ºQuartil: 118.6 97 RMax: Radia¸c˜ao Global M´axima Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 0.0 •M´edia amostral: 203.31 •M´aximo amostral: 988.0 •1ºQuartil: 0.1 •2ºQuartil: 14.1 •3ºQuartil: 347 VMed: Velocidade de Vento M´edia Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 0.0 •M´edia amostral: 2.016 •M´aximo amostral: 8.4 •1ºQuartil: 0.7 •2ºQuartil: 1.6 •3ºQuartil: 2.9 VMin: Velocidade de Vento M´ınima Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 0.0 •M´edia amostral: 0.3976 •M´aximo amostral: 3.3 •1ºQuartil: 0.0.0 •2ºQuartil: 0.3 •3ºQuartil: 0.6 98 VMax: Velocidade de Vento M´axima Descri¸c˜ao Num´erica Descri¸c˜ao Gr´afica •M´ınimo amostral: 0.0 •M´edia amostral: 5.35 •M´aximo amostral: 25.6 •1ºQuartil: 2.1 •2ºQuartil: 4.2 •3ºQuartil: 7.8 Seguidamente s˜ao consideradas duas medidas de inferˆencia estat´ıstica ´uteis para fazer considera¸c˜oes te´oricas sobre a vari´avel que est´a a ser estudada (seja Xessa vari´avel, onde X=x1, x2, ..., xn´e a amostra analisada): •Coeficiente de Curtose dado por γ2=µ4 µ2 2−3 = E[(X−E[X])4] (E[(X−E[X])2])2−3, onde µ2eµ4s˜ao, respetivamente, os momentos centrados de segunda e quarta ordens. Caso γ2>0, a distribui¸c˜ao da vari´avel ´e apelidada de leptoc´urtica; caso γ2<0, diz-se platic´urtica e apresenta uma distrubui¸c˜ao plana dos dados; quando γ2= 0, diz-se mesoc´urtica e deve ter um comportamento semelhante a uma popula¸c˜ao distribu´ıda segundo uma Normal. •Coeficiente de Assimetria dado por γ1=µ3 µ( 23 2)=E[(X−E[X])3] (E[(X−E[X])2])3 2 . γ1>0 significa que a popula¸c˜ao em an´alise apresenta uma m´edia superior `a mediana, denominando-se como assim´etrica `a direita. Caso contr´ario, ´e assim´etrica `a esquerda e apresenta maior mediana do que m´edia. µ3´e o momento centrado de terceira ordem. Vari´avel TMed TMin TMax HMed HMin HMax RMed γ10.054 0.056 0.055 −0.449 −0.382 −0.499 1.636 γ2−0.821 −0.819 −0.812 −0.219 −0.240 −0.191 1.542 Vari´avel RMin RMax VMed VMin VMax γ12.224 1.423 1.086 1.684 0.987 γ24.38 0.868 0.677 3.323 0.359