Previsão probabilística dos preços do mercado de eletricidade

Faculdade de Engenharia da Universidade do Porto Previsão Probabilística dos Preços do Mercado de Eletricidade Filipe António Sobral Sacramento Oliveira VERSÃO FINAL Dissertação realizada no âmbito do Mestrado Integrado em Engenharia Electrotécnica e de Computadores Major Energia Orientador: Ricardo Jorge Gomes Sousa Bento Bessa (Dr.) Co-orientador: José Nuno Moura Marques Fidalgo (Professor, Dr.) 16 de Fevereiro de 2015 © Filipe António Sobral Sacramento Oliveira, 2015 iii v Resumo Este documento tem por objetivo apresentar o trabalho realizado para obter modelos de previsão probabilística, bem como as respetivas previsões de preços de energia elétrica no Mercado Ibérico de Eletricidade (MIBEL). Os modelos de previsão são obtidos utilizando algoritmos baseados em Gradient Boosting. Estes algoritmos para além de previsões permitem obter estimativas das margens de incerteza. O Gradient Boosting é baseado no princípio de que se pode construir um estimador forte a partir da junção de vários estimadores fracos. As variáveis estudadas incluem valores históricos de preços, previsões meteorológicas, previsão de carga, entre outras. No final são avaliados os resultados de cada modelo para aferir quais as variáveis que afetam o preço, e qual o efeito que estas têm nas previsões. Este trabalho foi desenvolvido em ambiente R, e portanto foi necessário implementar programas em R para tratamento da informação e para automatização dos testes de previsão. A qualidade dos resultados obtidos permitiu confirmar o potencial do Gradient Boosting na previsão probabilística dos preços da energia elétrica. vii Abstract This document aims to present the work done in order to obtain forecasting models, as well as the respective forecasts of electricity prices in the Iberian Electricity Market (MIBEL). The forecast models are obtained using algorithms based on Gradient Boosting. These algorithms allow building models that can predict values from observations of selected explanatory variables. The variables studied include historical price values, weather forecasts, load forecasting, among others. At the end, results of each model are evaluated in order to assess which variables affect the price, and what effect they have on forecasts. One goal of this work was to learn and use statistical computing software R, so it was included in this document an attachment with information regarding the functions available in R for building models, to calculate forecasts and to create the graphs illustrating this document. ix Agradecimentos Quero em primeiro lugar agradecer à minha família, em especial à minha mãe que sempre me apoiou incondicionalmente em todos os momentos e em todas as decisões. Expresso também um sincero agradecimento ao Professor Doutor José Nuno Fidalgo e ao Doutor Ricardo Jorge Bessa pela disponibilidade demonstrada e pelo conhecimento e experiência partilhados. Um grande Obrigado a todos os meus amigos, da FEUP e de há muitos anos. Pela amizade, pelo apoio e pelos momentos passados. xvii Lista de tabelas Tabela 3.1 - Descrição das variáveis selecionadas para cada modelo ........................ 19 Tabela 3.2 – Variáveis selecionadas pelo gradient boost para o modelo 5 ................... 20 Tabela 4.1 – Correlações entre preço e previsões meteorológicas ............................ 26 Tabela 4.2 – Correlações entre preço e potência eólica produzida e previsão de carga ... 26 Tabela 4.3 – Descrição das variáveis selecionadas para cada modelo ........................ 27 Tabela 4.4 – Variáveis selecionadas pelo GLM Boost para o modelo 4 ........................ 33 xix Abreviaturas e Símbolos Lista de abreviaturas AR Autorregressiva ARIMA AutoRegressive Integrated Moving Average ARMA AutoRegressive Moving Average CRPS Continuous Rank Probability Score FEUP Faculdade de Engenharia da Universidade do Porto GAM Generalized Additive Model GD Gradiente Descendente GLM Generalized Linear Model MA Moving Average SEE Sistema Elétrico de Energia PRE Produção em Regime Especial PRO Produção em Regime Ordinário Lista de símbolos 𝑔(.) Estimador 𝜌 Função custo 𝑈 Vetor Gradiente 𝜐 Largura de passo de iteração 1 Introdução 1.1 Motivação Devido ao papel que os mercados de eletricidade têm vindo a assumir no panorama internacional, com cada vez maior relevância, surgiram zonas de mercado com países a comprar e a vender energia entre si, como são os casos do NordPool na Escandinávia, o BETTA (British Electricity Trading and Transmission Arrangements) e do MIBEL (Mercado Ibérico de Eletricidade). No entanto e apesar desta estabilidade global dos preços, estes mercados operam num regime de leilão com ofertas de compra e venda, e portanto torna-se necessário para os agentes que apresentam estas propostas, ajustarem-nas em função do preço espectável de mercado para que estas propostas sejam aceites e lhes possibilitem rendimentos otimizados. Assim, é necessário recorrer a modelos capazes de fornecer estimativas precisas dos preços da energia para cada período em negociação. Algoritmos de previsão probabilística têm vindo a ser utilizados na previsão de apostas desportivas, previsão meteorológica, previsão de produção de energia ou até de consumos energéticos. Os bons resultados obtidos noutras áreas sugerem que as técnicas utilizadas possam ser adaptadas ao caso dos mercados de eletricidade. Estes algoritmos também possibilitam a determinação de intervalos de confiança para as previsões, o que constitui uma vantagem relevante em diversas operações de mercado, porque permite obter uma medida do risco envolvido. É especialmente importante ser capaz de prever o comportamento dos preços do MIBEL até porque está anunciado o reforço das ligações elétricas entre Espanha e França, e esta medida implicará novas condições de exploração no MIBEL, permitindo assim a exportação de mais eletricidade de fontes renováveis para o resto da Europa e a importação de eletricidade mais barata. 2 Introdução 1.2 Objetivos Com este trabalho pretende-se criar modelos capazes de calcular previsões probabilísticas de preços do MIBEL a curto prazo, com a precisão necessária para que estas possam servir de apoio à decisão no processo de apresentação de propostas ao operador de mercado. Pretendese, além disso, estimar os intervalos de confiança associados às curvas de preços previstas. Os preços de mercado dependem de diversas condicionantes, como a disponibilidade e o preço de combustíveis, da potência instalada e da produção efetiva de eletricidade a partir de fontes renováveis de fontes renováveis, do diagrama de cargas das regiões em que o mercado opera, do histórico de preços de fecho de mercado, de acontecimentos esporádicos (tais como manutenções ou saídas de serviço forçadas), das estratégias de negócio dos agentes envolvidos, etc. Para a criação dos modelos é fundamental identificar as variáveis que potencialmente mais influenciam o preço de fecho de mercado, para que haja um maior cuidado no tratamento destas, de modo a aproveitar toda a sua capacidade discriminatória. Por outro lado, interessa identificar as variáveis que possam ser desprezadas pois apenas aumentam a incerteza do modelo para além de o tornarem mais lento. 1.3 Estrutura do trabalho Este documento encontra-se dividido em 5 capítulos mais um anexo com informações relativas à utilização da linguagem R. O Capítulo 1 é dedicado à introdução ao trabalho e ao tema em estudo, subdividindo-se este em Motivação, Objetivos e Estrutura do Trabalho. O Capítulo 2 é dedicado ao Estado da Arte. Na introdução a este capítulo é apresentada a perspetiva atual da utilização de previsões probabilísticas de preços e de outras variáveis. Seguidamente são apresentados os modelos de previsão encontrados na bibliografia consultada. No Capítulo 3 é apresentada a metodologia utilizada durante este trabalho. A teoria que serviu de base a todo o trabalho realizado encontra-se descrita neste capítulo. No final encontra-se o resumo dos passos seguidos para se obterem resultados e para que estes fossem avaliados. O Capítulo 4 é onde se encontram os resultados do estudo preliminar realizado com os dados iniciais, os resultados das previsões calculadas e os resultados da avaliação destas. No final deste capítulo é feita a comparação entre os modelos obtidos. As conclusões e trabalho futuro encontram-se no Capítulo 5. Neste capítulo são tiradas conclusões sobre os objetivos do trabalho, a qualidade das previsões calculadas e são dadas sugestões sobre como melhorar os resultados obtidos. 3 1.3 Estrutura do trabalho A última parte deste documento é constituída pelas referências bibliográficas consultadas e pelos Anexos. Neste último encontra-se a descrição do que foi feito em linguagem R, ilustrado com exemplos de código implementado e comentários e explicações sobre a linguagem de programação, o ambiente de programação e sobre as funções utilizadas. 5 Estado da Arte Neste capítulo serão expostos os métodos atualmente utilizados para previsão dos preços de eletricidade. 2.1 Introdução Em termos históricos a previsão de preços de eletricidade começou por ser feita com recurso a técnicas regressivas, estas na sua forma clássica são baseadas no critério dos mínimos quadrados no qual o modelo é calculado por forma a que a soma das diferenças entre os valores observados e os valores previstos seja minimizada. Mesmo hoje em dia, e apesar do grande número de alternativas que já existem, os métodos de regressão linear continuam a ser muito utilizados para previsões de preços de energia. A partir da década de 60 com a evolução da complexidade dos SEE foi necessário recorrer a metodologias de previsão cada vez mais elaboradas. O desenvolvimento do computador pessoal a partir da década de 80 permitiu a exploração de algoritmos de previsão mais complexos como sejam as Redes Neuronais Artificiais, Fuzzy Logic e Algoritmos Genéticos. No artigo [1] de 2011 apesar de não se tratar a problemática da previsão de preços de eletricidade, são calculadas previsões de produção de energia eólica recorrendo a uma metodologia muito semelhante à utilizada neste trabalho para fazer previsões de preços. A aplicação de previsão probabilística não surge frequentemente na literatura referente a previsão de preços da eletricidade. Só recentemente começam a surgir trabalhos como [2] e [3]. Nestes artigos os resultados obtidos revelam-se bastante promissores quando comparados com resultados obtidos por outros métodos, tendo sido registadas previsões mais precisas com regressões por quantis do que com outros métodos regressivos não paramétricos. 12 Metodologia estações meteorológicas e ainda a previsão de produção eólica total estão correlacionadas com o preço da eletricidade, foi calculada a correlação de Spearman de cada uma destas variáveis com o preço. Os resultados das correlações foram incluídos no capítulo 4, onde são apresentados os resultados deste trabalho. 3.1.1. Tipo de previsão A previsão de preços é realizada sob a forma de intervalos onde o preço verificado se encontrará com uma determinada probabilidade. Os intervalos considerados são quantis de 5% até 95%, com uma amplitude de 5%, o que resulta em 19 intervalos. O Gráfico 3.1 representa um exemplo do resultado da previsão de preços. Neste a linha a preto representa o preço real que se verificou no instante de tempo representado, e em tons de cinzento estão representadas as previsões. Gráfico 3.1 – Exemplo de previsão de preço 3.2 Técnicas de regressão 3.2.1. Gradient Boosting O Boosting é um algoritmo de aprendizagem automática caracterizado por iterativamente fazer a aprendizagem de estimadores fracos e adicioná-los por forma a obter um estimador forte. Um estimador fraco é definido pela sua fraca correlação com o que se pretende classificar, o que significa que o não é possível fazer previsões com rigor recorrendo apenas a 13 3.2 Técnicas de regressão este tipo de variável explicativa. Pelo contrário um estimador forte é aquele que está bem correlacionado com o que se pretende prever. O gradient boosting é uma técnica que permite a construção de modelos preditivos baseados em dados fornecidos como entradas. Esta técnica é tipicamente utilizada em problemas que envolvam regressões. Especificamente neste trabalho é utilizado o gradient boosting com uma função de custo relacionada com a regressão por quantis, que será descrita em pormenor mais adiante neste capítulo. Este método foi inventado por Jerome H. Friedman e descrito em dois artigos publicados em fevereiro e março de 1999, ver [21] e [22], sendo que o primeiro descreve o algoritmo e o segundo artigo descreve como melhorar a performance do método. O método do gradient boosting baseia-se na minimização de uma função de custo que penaliza a diferença entre os valores obtidos pelo modelo preditivo e os valores medidos. Este processo é aplicado para cada variável explicativa selecionada obtendo-se um estimador por cada uma destas variáveis. Após este processo os estimadores são combinados resultando numa função de estimadores, a qual será representativa do modelo. De seguida descreve-se um procedimento generalista de como o processo se desenrola para a obtenção desta função. Partindo de um conjunto de dados (X1, Y1), …, (Xn,Yn) escolhe-se um procedimento de base, que no caso do glmboost é um modelo linear e no caso do gamboost é um spline cúbico, para obter um estimador 𝑔. (𝑋1,𝑌1)𝑃𝑟𝑜𝑐𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 𝑏𝑎𝑠𝑒 → 𝑔(.) O processo é repetido para as outras variáveis relevantes para o problema: 𝐷𝑎𝑑𝑜𝑠 𝑃𝑜𝑛𝑑𝑒𝑟𝑎𝑑𝑜𝑠 1𝑃𝑟𝑜𝑐𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 𝑏𝑎𝑠𝑒 → 𝑔[1](.) 𝐷𝑎𝑑𝑜𝑠 𝑃𝑜𝑛𝑑𝑒𝑟𝑎𝑑𝑜𝑠 2𝑃𝑟𝑜𝑐𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 𝑏𝑎𝑠𝑒 → 𝑔[2](.) … … 𝐷𝑎𝑑𝑜𝑠 𝑃𝑜𝑛𝑑𝑒𝑟𝑎𝑑𝑜𝑠 𝑀𝑃𝑟𝑜𝑐𝑒𝑑𝑖𝑚𝑒𝑛𝑡𝑜 𝑏𝑎𝑠𝑒 → 𝑔[𝑀](.) No final os estimadores são agregados numa única expressão: 𝑓(.)= ∑𝛼𝑚∗𝑔[𝑀] 𝑀 𝑚=1 (3.1) O termo “dados ponderados” significa que são atribuídos pesos individuais (𝛼) a cada uma das variáveis. De seguida é apresentado o algoritmo do gradiente descendente funcional, ou gradient boost, cujo nome deriva da demonstração de Breiman, ver [23] e [24] de que o método original 14 Metodologia AdaBoost poderia ser representado por um algoritmo de “descida mais íngreme”, Steepest Descent. A estimação de 𝑓∗(.) com boosting pode ser feita considerando a função risco empírico dada por 𝑛−1 ∗∑𝜌(𝑌𝑖,𝑓(𝑋𝑖)) 𝑛 𝑖=1 que representa a média da função custo, seguindo iterativamente pela descida mais íngreme dentro do espaço da função. Assim Friedman, ver [21], propôs o seguinte algoritmo: 1. Inicializar 𝑓[0] com um valor inicial, escolhas mais usuais podem ser: a. 𝑓[0] ≡argmin𝑛−1∑𝜌(𝑌𝑖,𝑐) 𝑛 𝑖=1 b. 𝑓[0] ≡0 Como esta é a iteração 0, definir 𝑚=0. 2. Especificar um c1onjunto de estimadores de base. Estimadores de base são estimadores de regressão simples com um conjunto de variáveis de entrada simples e resposta univariada. Os conjuntos de variáveis de entrada podem ser diferentes entre estimadores de base.1 Normalmente as variáveis de entrada dos estimadores de base são pequenos subconjuntos do grupo das variáveis preditoras 𝑥1,… 𝑥1𝑛. 3. Incrementar m de 1. Calcular o vetor gradiente negativo −𝜕 𝜕𝑓𝜌(𝑌,𝑓) e avaliar em 𝑓[𝑚−1](𝑋𝑖): 𝑈𝑖=− 𝜕 𝜕𝑓𝜌(𝑌𝑖,𝑓)|𝑓=𝑓[𝑚−1](𝑋𝑖),𝑖 =1,…,𝑛 (3.2) 4. Aplicar cada um dos estimadores de base ao vetor gradiente negativo, isto é, usar cada um dos estimadores da regressão especificados no passo 2 separadamente no vetor gradiente negativo. 5. Selecionar o estimador de base que melhor se adequa a 𝑈[𝑚] de acordo com o critério da soma residual dos quadrados e estabelecer Û[𝑚] igual aos valores do estimador de base mais adequado 6. Atualizar 𝑓[𝑚](.)=𝑓[𝑚−1](.)+𝜐∗𝑔[𝑚](.), onde 0<𝜐≤1 é um fator de comprimento de passo de iteração 7. Iterar os passos 2 a 6 até se verificar 𝑚=𝑚𝑠𝑡𝑜𝑝 sendo 𝑚𝑠𝑡𝑜𝑝 a iteração de paragem. O principal parâmetro de afinação do método é a iteração de paragem, a qual pode ser determinada por validação cruzada. Mais à frente neste documento será exposto o procedimento que levou à determinação deste parâmetro. A importância da escolha deste parâmetro reside no seguinte facto: se a iteração de paragem for demasiado pequena, correse o risco do algoritmo introduzir poucas variáveis no modelo; se for demasiado grande, para além do problema óbvio de maior esforço computacional, acresce um problema maior que é a ocorrência de overfitting, que significa que o modelo se tornou demasiado adaptado aos dados 15 3.2 Técnicas de regressão do conjunto de treino e portanto não é capaz de extrapolar resultados coerentes quando lhe são apresentados na entrada os dados do conjunto de validação. A escolha da largura do passo, 𝜐 no passo 4 do algoritmo não é muito importante desde que este seja pequeno, como 𝜐=0.1 por exemplo, [25], [26] e [27]. Por outro lado, um valor mais pequeno de 𝜐 tipicamente implica um maior número de iterações, o que por sua vez leva a um aumento do tempo de cálculo. A validação cruzada permite avaliar a capacidade de generalização de um modelo a partir de um conjunto de dados. O conceito geral é a partição do conjunto de dados em subconjuntos mutuamente exclusivos, e posteriormente utilizar alguns destes subconjuntos para a estimação dos parâmetros do modelo, e os restantes subconjuntos para a validação do modelo. Para este trabalho utilizou-se o método k-fold que consiste em dividir os dados em k subconjuntos sendo um destes para validação, e os restantes k-1 para estimação dos parâmetros do modelo. Este processo é realizado k vezes alternando ciclicamente o subconjunto de teste. No final calculase a precisão do modelo pela minimização da função custo. Na realização deste processo os dados foram divididos em 5 subconjuntos, e portanto k=5. Constatou-se que a partir das 5000 iterações não ocorriam melhorias significativas da precisão do modelo, e portanto ficou definido que a iteração de paragem seria 𝑚𝑠𝑡𝑜𝑝 =5000. Nos gráficos Gráfico 3.2 e Gráfico 3.3 estão representados os resultados da validação cruzada para 5000 e 15000 iterações para o conjunto de teste. Gráfico 3.2 - Validação Cruzada com 5000 iterações Gráfico 3.3 - Validação Cruzada com 15000 iterações Analisando os gráficos é possível verificar que até às 15000 iterações não existe overfitting, pois se houvesse o gráfico seria crescente a partir da iteração onde este começasse indicando 16 Metodologia que a função de custo (indicador de desempenho) seria crescente e portanto a capacidade de generalização do modelo estaria a diminuir. 3.2.2. Regressão por quantis A regressão por quantis é utilizada no gradient boosting por via da função de custo e da sua derivada para criar a previsão probabilística pretendida com este trabalho. Existem aplicações em que, para além de ser necessário estimar o valor de uma grandeza, é também necessário estimar a incerteza que lhe está associada. A análise da distribuição do erro permite estimar a incerteza, mas para efetuar essa análise é necessário partir de pressupostos sobre a distribuição que o erro adota. No entanto existem tipos de previsões que não se encaixam devidamente nos tipos de distribuição assumidos. Assim pode ser vantajoso tratar a incerteza do modelo de forma separada e sem assumir algum tipo de distribuição logo à partida, e portanto não influenciar nem os resultados das previsões nem o erro destas durante o processo de previsão. Para tal pode-se fazer uma regressão por quantis. Os quantis são uma medida estatística que quantifica um conjunto de dados. As linhas dos quantis separam os dados de modo a que o número de observações por baixo da linha corresponda a um determinado rácio que define o quantil. Existem quantis específicos como os tercis, que dividem os dados em 3 áreas cada uma contendo uma parte igual do total das observações, os quartis, que dividem as observações em 4 áreas, os quintis, os decis, os percentis, entre outros. Os quantis de uma amostra ou população são por vezes expressos na forma 𝜏∈[0,1] e portanto para o 𝜏-ésimo quantil, 100×𝜏% das observações deverão ter um valor inferior ao quantil 𝜏. Para encontrar o melhor valor que divide as observações é introduzida uma função de custo linear e assimétrica que penaliza os desvios de acordo com o lado da linha em que os dados são colocados. Assim sendo, minimizar o custo total é o mesmo que minimizar o erro total para um determinado quantil. Seja 𝑌 uma variável aleatória com função distribuição de probabilidade 𝐹𝑌, e seja 𝜏 um número real 0<𝜏≤1. Então o 𝜏-ésimo quantil de 𝐹𝑌 será: 𝑞𝑌(𝜏)=𝐹𝑌 −1(𝜏)=inf{𝑦:𝐹𝑌(𝑦)≥𝜏} (3.3) Isto significa que 100×𝜏% do valor de 𝑌 é inferior a 𝑞𝑌(𝜏). Isto quer dizer que o quantil 𝜏 para 𝑌 é uma função para a qual: 𝐹(𝑞𝑌(𝜏))=ℙ(𝑌<𝑞𝑌(𝜏))=𝜏 (3.4) Agora assumindo que existe um conjunto de dados (𝑌,𝑋) com N observações em que 𝑌 é 𝑁×1 e 𝑋 é 𝑁×𝑘, então temos o modelo linear: 𝑦𝑡=𝑄 (𝜏,𝑥𝑡)+𝜀𝑡=𝑥𝑡𝛽+𝜀𝑡 𝑡=1,…,𝑁 (3.5) 17 3.2 Técnicas de regressão A função de custo é dada por: 𝜌𝜏(𝜀)={(𝜏−1)𝜀, 𝑠𝑒 𝜀<0 𝜏𝜀, 𝑠𝑒 𝜀≥0 (3.6) E portanto o melhor 𝜏-ésimo estimador de 𝛽 pode ser determinado minimizando a função objetivo: 𝑉𝑁(𝛽;𝜏)=1 𝑁∑𝜌𝜏(𝑦𝑡−𝑥𝑡𝛽) 𝑁 𝑡=1 (3.7) Em problemas reais por vezes surgem processos que não são aproximáveis por relações lineares, e portanto é importante admitir algumas não linearidades no modelo. O gradient boosting é frequentemente utilizado com splines para modelizar estas não linearidades. Os splines são polinómios de grau 𝑚 em intervalos definidos por uma sequência de nós e que são 𝑚−1 vezes diferenciáveis nesses nós. Por cada intervalo entre dois nós existe um polinómio de grau m que minimiza o caminho entre esses nós. Existem duas categorias de splines, cada uma com o seu objetivo. Splines de interpolação que passam por todos os nós e splines de aproximação que passam perto dos nós. 3.3 GLM (Generalized Linear Model) Boosting Este método do Gradient Boosting tenta enquadrar os dados que lhe são fornecidos num modelo linear da forma: 𝑔(𝜇)=𝛽0+𝛽1𝑥1+⋯+𝛽𝑝𝑥𝑝 (3.8) Sendo que as variáveis explicativas do problema são da forma 𝑥=(𝑥1,…,𝑥𝑝) e os parâmetros 𝛽 representam os pesos associados a cada variável do modelo. Este tipo de modelos tende a produzir resultados pouco precisos devido ao facto de que por vezes a relação entre as variáveis explicativas e o que se pretende prever pode não ser linear. Apesar do exposto serão realizadas previsões utilizando este modelo para aferir a qualidade dos resultados obtidos e comparar com os resultados de outros modelos. 3.4 GAM (Generalized Aditive Model) Boosting Neste método ao contrário do anterior já não se assume à partida uma relação linear entre variáveis e grandeza a prever. Assim sendo, neste método são calculadas funções dependentes das variáveis que são agregadas na função descritiva do modelo em uso da seguinte forma: 𝑔(𝜇)=𝛽0+𝑓1+⋯+𝑓𝑝 (3.9) 18 Metodologia Onde 𝑓1, …,𝑓𝑝 são funções de 𝑥=(𝑥1,…,𝑥𝑝) respetivamente. As funções 𝑓𝑝 incluem funções lineares simples, bem como funções não lineares diferenciáveis, ou seja splines. É possível assim definir cada tipo de variável como linear, cíclica ou outra para melhor as adequar ao modelo. Para além dos efeitos que se podem atribuir às variáveis, pode também incluir-se a função 𝑔(𝜇) numa família de funções. A escolha de uma família tem consequências na função custo e no correspondente gradiente negativo, já que famílias diferentes têm funções custo, bem como gradientes diferentes. Uma escolha criteriosa da função custo permite a estimação de qualquer característica desejada da distribuição. Isto combinado com um grande número de estimadores garante um vasto conjunto de modelos que podem ser analisados com boosting. Por defeito a família utilizada é a gaussiana (Gaussian), mas neste caso foi escolhida a família de funções de regressão por quantis (QuantReg), que permite criar uma previsão probabilística. 3.5 Modelos Entradas/Saídas Para a obtenção dos modelos foi criado em R um programa que faz a importação dos dados dos ficheiros de texto onde estes se encontram, e recorrendo às funções glmboost ou gamboost (ver anexo) gera um modelo composto pelas variáveis explicativas definidas pelo utilizador. No entanto, em alguns modelos, nem todas as variáveis selecionadas à partida fazem parte do modelo. Isto acontece porque o gradient boost também tem um processo de seleção de variáveis, no qual de entre as variáveis que lhe são passadas pelo utilizador são escolhidas as variáveis que melhor descrevem o modelo. Nesta secção serão descritos os modelos testados neste trabalho. Foram criados 8 modelos diferentes, 4 com GLM e 4 com GAM. A diferença entre modelos reside não só no método de boosting na sua génese, mas também nas variáveis selecionadas para fazerem parte de cada modelo. A Tabela 3.1 dá uma visão global dos modelos. 19 3.5 Modelos Entradas/Saídas Variáveis explicativas Previsões Valores passados das séries Calendarização Carga, Eólica e Meteorológica Semana Anterior Dia Anterior 2 Horas Anteriores Mês Dia do mês Dia da Semana Hora do dia Modelos GLM Boost 1         2         3         4         GAMBoost 5         6         7         8         Tabela 3.1 - Descrição das variáveis selecionadas para cada modelo A variável “Semana Anterior” assume o valor do preço da energia 168 horas antes da previsão que se pretende calcular e a variável “Dia Anterior” assume o valor do preço da eletricidade 24 horas antes da previsão que se pretende. O modelo 1 é criado por GLM e as variáveis selecionadas foram as previsões meteorológicas para as 4 regiões consideradas, a previsão de carga para o dia seguinte e a previsão de produção de energia eólica. Para os modelos 2 e 6 foram selecionadas apenas as variáveis preço da eletricidade na semana anterior e preço da energia no dia anterior. Com estes modelos pretende-se avaliar a influência dos valores passados das séries de preços na previsão de preços futuros. A diferença entre estes modelos reside no método que lhes dá origem, o modelo 2 é criado por GLM e o modelo 6 é criado por GAM. Os modelos 3 e 7 são semelhantes aos modelos 2 e 6, mas foram acrescentadas as variáveis preço da hora anterior e os preços nas duas horas anteriores à primeira hora a prever. Estes modelos foram criados para perceber a influência da tendência de preços ao longo do horizonte de previsão pretendido para este trabalho. Estas variáveis permitem que o modelo para a previsão da hora t disponha da tendência que os preços vinham a seguir nas horas t-1 e t-2. No modelo 4 foram selecionadas as variáveis correspondentes a valores passados das séries de preços e as previsões do modelo 1 para determinar se a qualidade das previsões de preços melhora com a introdução de mais variáveis no modelo ou se, pelo contrário, o excesso de variáveis só provoca ruido nas previsões calculadas, aumentando assim o desvio destas em relação aos preços verificados. O modelo 5 foi criado para tirar proveito das capacidades do GAM boost na modelização de relações não lineares entre as variáveis explicativas e a grandeza a prever. Assim foram 20 Metodologia adicionadas 4 novas variáveis a este modelo comparativamente ao modelo 1. As variáveis são o mês do ano, o dia do mês, o dia da semana, e a hora do dia. O modelo 8 combina o pretendido com os modelos 4 e 5. No modelo 8 pretende-se tirar partido das variáveis mês, dia do mês, dia da semana e hora do dia, que não podiam entrar no modelo 4, para determinar se com mais variáveis é possível melhorar a qualidade das previsões obtidas. Mas como já foi referido nem todas as variáveis explicativas passadas ao GLM e ao GAM são selecionadas por estes métodos para fazer parte do modelo. Na Tabela 3.2 encontram-se as variáveis selecionadas pelo gradient boosting em diferentes quantis do modelo 5. Quantil 5% 25% 50% 75% 95% Previsão de Carga      Previsão de Produção Eólica      Previsão de Precipitação em Alto Douro Previsão de velocidade de vento em Alto Douro     Previsão de Temperatura em Alto Douro    Previsão de Precipitação em Ameixeiras Previsão de velocidade de vento em Ameixeiras     Previsão de Temperatura em Ameixeiras  Previsão de Precipitação em Lourinhã Previsão de velocidade de vento em Lourinhã Previsão de Temperatura em Lourinhã  Previsão de Precipitação em Torre Miró Previsão de velocidade de vento em Torre Miró     Previsão de Temperatura em Torre Miró  Mês     Dia do mês     Dia da semana Hora     Tabela 3.2 – Variáveis selecionadas pelo gradient boost para o modelo 5 Como se pode verificar pela análise da tabela, a previsão de precipitação não é selecionada para fazer parte do modelo 5. O GAM boost seleciona todas as previsões de temperatura, embora em quantis diferentes, o que indica que esta variável é mais explicativa do preço do que a previsão de precipitação. A produção de eletricidade em aproveitamentos hidroelétricos influencia o preço de mercado da eletricidade. No entanto a maioria da água turbinada nestes aproveitamentos pode demorar vários dias desde que se precipita sob a forma de chuva até chegar às albufeiras das barragens. Neste caso apenas é considerada a previsão de precipitação 21 3.7 Resumo para o dia seguinte, pelo que um período de 24 horas no máximo é pouco tempo para que a precipitação que se prevê cair no terreno chegue às albufeiras. Além disso, a produção hídrica é usada frequentemente para compensação ou reserva, pelo que o facto de existir água disponível não implica que esta venha a ser turbinada no imediato. Nesta tabela observa-se que o dia da semana não é selecionado, o que não era espectável, dado ser uma variável relevante em termos de previsão de carga. Para tentar esclarecer esta questão, experimentou-se calcular os modelos mas sem a variável Previsão de Carga e verificou-se que, nesta situação, a variável Dia da Semana passa a ser selecionada. Conclui-se então que o Dia da Semana não foi selecionado anteriormente porque os efeitos desta variável se encontram de certo modo representados na variável Previsão de Carga. Ou seja, o Dia da Semana efetivamente influencia a carga e que esta, por sua vez, influencia o preço. O programa criado também calcula as previsões, para tal utiliza a função predict do R. Esta função recebe como argumentos o modelo com o qual se pretendem fazer as previsões e o conjunto de dados de validação. Os valores previstos são organizados numa matriz em que cada coluna corresponde a um quantil e cada linha corresponde à previsão de uma hora. Cada modelo deu origem a uma matriz com 19 colunas, correspondentes aos quantis desde 5% até 95% em intervalos de 5%, e 6000 linhas ou previsões hora a hora. No entanto é de salientar que neste trabalho são feitas previsões para um horizonte máximo de 24 horas, pelo que as linhas que compõem as matrizes referidas estão organizadas em 150 conjuntos de 24 previsões, cada um correspondendo a um dia do conjunto de dados de validação. 3.6 Avaliação dos Resultados Para avaliar os resultados obtidos foram usadas três métricas, sendo estas a calibração dos dados, a largura dos intervalos e CRPS. Estas métricas foram calculadas recorrendo a três funções do R que entre outros, recebiam como parâmetros as previsões realizadas pelos modelos testados e os valores dos preços que se verificaram no MIBEL no período para o qual se fizeram as previsões, que seriam usados como referência. No final foram feitos gráficos comparativos dos resultados das avaliações dos diferentes modelos. Um requisito para previsões probabilísticas é que as probabilidades nominais (ou proporções nominais) das previsões dos quantis correspondam às probabilidades efetivas. Posto de outra forma, a calibração avalia a diferença entre os quantis previstos e os quantis nominais que representam. A amplitude dos intervalos é a tendência das previsões probabilísticas para previsões discretas medida pelo tamanho médio dos intervalos de previsão (distância entre quantis). Os quantis são reunidos aos pares para se obterem intervalos com diferentes taxas de cobertura. Isto dá uma indicação do nível de utilidade onde são necessários intervalos estreitos. Esta métrica não depende das observações. 28 Resultados Obtidos Gráfico 4.1 – Previsão com o modelo 6 para 28/4/2014 No Gráfico 4.2 estão representadas as previsões obtidas com o modelo 8 para o dia 28 de abril de 2014. Gráfico 4.2 - Previsão com o modelo 8 para 28/4/2014 Mais adiante neste capítulo serão avaliados os resultados de cada modelo em pormenor, mas numa primeira análise é possível concluir que em geral o modelo 8 gera previsões em 29 4.2 Previsões quantis com menor amplitude do que o modelo 6. Devido a este facto é mais frequente que os quantis das previsões do modelo 8 não contenham o preço da eletricidade do que os quantis do modelo 6. Contudo não se pode concluir que o modelo 6 produza previsões com mais qualidade do que o modelo 8 apenas pela análise visual destes gráficos. 4.3 Avaliação dos Resultados A primeira avaliação realizada às previsões obtidas foi a análise do erro médio absoluto. No Gráfico 4.3 encontra-se representada a evolução ao longo do horizonte de previsão do erro médio absoluto para o quantil 50% do modelo 8 para as previsões calculadas desde janeiro de 2014 até setembro de 2014. Gráfico 4.3 - Erro médio absoluto das previsões hora a hora para o quantil 50% do modelo 8 Uma razão para a tendência crescente do erro com o horizonte de previsão será descrita mais adiante neste capítulo. Por inspeção visual do gráfico pode-se constatar que o máximo que o erro médio atinge é de aproximadamente 7€/MWh para este modelo. O modelo 7 tem o máximo do erro médio absoluto do quantil 50% para aproximadamente 12€/MWh, e o modelo 5 para o mesmo quantil tem no máximo 9.7€/MWh de erro médio absoluto. O erro absoluto médio obtido para o modelo 8 é inferior ao registado em outros estudos similares, como o [4]. Embora o conjunto de teste utilizado em [4] não seja coincidente com o 30 Resultados Obtidos utilizado neste estudo, o desempenho obtido para o quantil 50% mostra que a técnica aplicada tem potencial, além de permitir estimar também a incerteza associada às previsões. 4.3.1. Resultados do GLM Boost Os resultados das previsões foram posteriormente avaliados segundo as três métricas já referidas no capítulo 3, a calibração, a amplitude dos intervalos e o CRPS. Com o resultado das avaliações foram construídos os gráficos comparativos Gráfico 4.4, Gráfico 4.5 e Gráfico 4.6. O Gráfico 4.4 apresenta o desvio de calibração obtido para os modelos construídos com GLM, a linha de referência representa a calibração perfeita. Gráfico 4.4 – Desvio de calibração dos modelos produzidos com GLM Boost Analisando o Gráfico 4.4 é possível constatar que o modelo 1, que é constituído pelas previsões meteorológicas pela previsão de carga e pela previsão de produção eólica, apresenta bons resultados para os quantis das extremidades. No entanto para quantis intermédios apresenta os piores resultados de entre os modelos testados. O modelo 2 apresenta uma boa calibração e é apenas constituído pelas variáveis semana anterior e dia anterior, pelo que necessita de pouco esforço computacional. Analisando apenas a calibração não seria previsível que um modelo que conta com tão poucas variáveis explicativas para o descrever apresentasse uma calibração tão próxima da calibração perfeita. No entanto a qualidade de um modelo deve ser avaliada de uma perspetiva mais alargada e não apenas pelo resultado de uma única métrica. Analisando o gráfico pode-se verificar que o modelo 3 apresenta uma melhor calibração até ao quantil 45%, a partir do qual é o modelo 4 o que apresenta a melhor calibração entre estes. 31 4.3 Avaliação dos Resultados Ao contrário do modelo 2, os modelos 3 e 4 apresentam resultados que não eram previsíveis, seria de esperar que no caso do modelo 3 com a introdução das variáveis preço da energia na hora anterior e há duas horas atrás os resultados fossem ainda melhores que os obtidos para o modelo 2. O mesmo seria de esperar para o modelo 4 que para além das variáveis do modelo 3 conta ainda com as previsões meteorológicas. É necessário portanto analisar os resultados das outras métricas. A largura dos intervalos indica a distância entre dois quantis diferentes. No Gráfico 4.5 está representada a amplitude dos intervalos. Gráfico 4.5 – Amplitude dos intervalos interquantis dos modelos produzidos com GLM Boost Analisando este gráfico torna-se evidente a razão pela qual o modelo 2 apresenta resultados tão bons para a calibração. A razão é que devido à amplitude dos intervalos de previsão ser mais elevada do que nos modelos 3 e 4, as previsões calculadas correspondem mais frequentemente aos seus respetivos quantis. Já nos modelos 3 e 4 ocorre a situação oposta, como os intervalos interquantis são menores é mais frequente a ocorrência de previsões que não cabem nestes intervalos piorando assim o resultado da calibração. O modelo 1 também apresenta uma amplitude muito elevada, quatro vezes superior à do modelo 4. O CRPS mede o desempenho geral da previsão probabilística. No Gráfico 4.6 estão representados os resultados do CRPS para os modelos de 1 a 4 em função do horizonte de previsão. 32 Resultados Obtidos Gráfico 4.6 – Resultados do CRPS dos modelos produzidos com GLM Boost Analisando o Gráfico 4.6 é possível verificar que efetivamente o modelo que produz melhores resultados é o 4. Apesar de este começar com bons resultados para as primeiras horas do dia, estes vão em geral piorando com o horizonte de previsão. Isto acontece pois para prever a primeira hora as variáveis hora anterior e duas horas antes são os preços verificados nas duas últimas horas do dia anterior, pelo que o modelo segue a tendência dessas duas horas anteriores. Nas previsões seguintes são atribuídos a estas variáveis valores de previsões passadas, pelo que ao longo do horizonte de previsão ocorre propagação do erro. A ocorrência de propagação do erro ocorre devido à utilização recursiva de previsões em instantes passados para calcular novas previsões, como as previsões calculadas têm desvios relativamente ao preço observado, ao calcular novas previsões recorrendo a dados com erros propaga-se o erro para a previsão seguinte. O modelo 3 também conta com as variáveis referentes às duas horas anteriores pelo que apresenta um comportamento semelhante ao do modelo 4, mas por não terem sido usadas as variáveis referentes às previsões o modelo 3 produz resultados piores comparativamente aos do modelo 4. O modelo 2 apesar dos bons resultados de calibração produz em geral piores resultados que o modelo 4, no entanto a partir da hora 18 produz os melhores resultados de entre os modelos comparados. Uma possível explicação para este comportamento é que o acumular de erros de previsão para a hora seguinte nos modelos 3 e 4 produz variáveis hora antes e duas horas atrás com valores muito diferentes dos valores verificados, o que introduz erros cada vez maiores na previsão da hora seguinte. Apesar do exposto o modelo 4 é aquele que produz os melhores resultados de entre os apresentados porque para além de integrar os 33 4.3 Avaliação dos Resultados valores passados das séries de preços, inclui também as previsões de carga, de produção eólica e meteorológicas. As variáveis que o GLM selecionou para o modelo 4 encontram-se resumidas na Tabela 4.4. Quantil 5% 25% 50% 75% 95% Previsão de Carga      Previsão de Produção Eólica      Previsão de velocidade de vento em Alto Douro      Previsão de Temperatura em Alto Douro      Previsão de velocidade de vento em Ameixeiras      Previsão de Temperatura em Ameixeiras      Previsão de velocidade de vento em Lourinhã      Previsão de Temperatura em Lourinhã      Previsão de Precipitação em Torre Miró      Previsão de velocidade de vento em Torre Miró      Previsão de Temperatura em Torre Miró      Preço na Semana Anterior      Preço no Dia Anterior      Preço na hora anterior      Preço há duas horas atrás      Tabela 4.4 – Variáveis selecionadas pelo GLM Boost para o modelo 4 4.3.2. Resultados do GAM Boost Os resultados da avaliação das previsões produzidas utilizando os modelos obtidos com GAM são em geral semelhantes aos resultados do GLM. No Gráfico 4.7 estão representados os resultados do desvio de calibração. 34 Resultados Obtidos Gráfico 4.7 – Desvio de calibração dos modelos produzidos com GAM Boost Analisando o Gráfico 4.7 verifica-se que as previsões obtidas têm uma calibração muito semelhante à calibração do GLM. O modelo 5, tal como o modelo 1 apresenta bons resultados para os quantis das extremidades mas tem os piores resultados de entre os modelos testados para valores intermédios dos quantis. Neste caso o modelo 6 também apresenta muito bons resultados para a calibração mas, como se poderá observar adiante neste documento, tal como o modelo 2 tem a amplitude dos intervalos elevada o que justifica a boa calibração. Os modelos 7 e 8 produzem resultados de calibração pouco diferentes dos resultados dos modelos 3 e 4. No Gráfico 4.7 pode-se observar ainda que os modelos que apresentam uma maior diferença entre GLM e GAM são os modelos 1 e 5. A diferença entre a calibração dos resultados do modelo 5 e a referência é superior à diferença entre a calibração dos resultados do modelo 1 e a referência. O desvio máximo da calibração do modelo 1 é 0,177 para o quantil 20% e o desvio máximo da calibração do modelo 5 é de 0,287 para o quantil 40%. A amplitude dos intervalos, tal como acontece para a calibração, é muito semelhante entre modelos produzidos por GAM e GLM. No Gráfico 4.8 está representada a amplitude dos intervalos dos modelos produzidos com GAM. 35 4.3 Avaliação dos Resultados Gráfico 4.8 - Amplitude dos intervalos interquantis dos modelos produzidos com GAM Boost Uma vez mais se pode verificar que no modelo constituído por semana e dia anteriores, agora para o modelo 6, a amplitude dos intervalos é grande quando comparada à dos modelos 7 e 8, este facto explica a razão dos bons resultados da calibração do modelo 6 face aos modelos 7 e 8 que são mais completos em termos de variáveis explicativas. Nos resultados do CRPS já é possível observar diferenças mais significativas. Os resultados do CRPS encontram-se no Gráfico 4.9. Gráfico 4.9 - Resultados do CRPS dos modelos produzidos com GLM Boost 36 Resultados Obtidos A primeira grande diferença é a melhoria dos resultados do modelo 5 face ao modelo 1. Para o modelo 5 para além das variáveis explicativas selecionadas para o modelo 1 selecionaram-se também as variáveis mês do ano, dia do mês, dia da semana e hora do dia. Estas variáveis têm influência no preço final pois por exemplo em horas de cheia do diagrama de cargas o preço da energia é mais elevado do que nos períodos de vazio, pelo que a varável hora do dia é de grande importância quando se pretende analisar o preço da energia ao longo do dia. Uma explicação análoga pode ser considerada para a variável mês já que principalmente no verão e no inverno o preço da eletricidade é afetado pelo peso no diagrama de cargas dos dispositivos de climatização. Os modelos 2 e 6 têm resultados muito semelhantes. Os modelos 3, 4, 7 e 8 vão ser objeto de uma comparação mais detalhada nos gráficos 4.10, 4.11 e 4.12. 37 4.3 Avaliação dos Resultados Gráfico 4.10 - Desvio de calibração dos modelos 3, 4, 7 e 8 Gráfico 4.11 - Amplitude dos intervalos interquantis dos modelos 3, 4, 7 e 8 Gráfico 4.12 - Resultados do CRPS dos modelos 3, 4, 7 e 8 Analisando os gráficos é possível verificar que a calibração é semelhante quer para os modelos obtidos por GAM quer para os modelos obtidos por GLM. No entanto os modelos obtidos por GLM têm a menor amplitude dos intervalos. Na apreciação global fornecida pelo CRPS podese verificar que comparando modelos equivalentes, modelos 3 e 7 e modelos 4 e 8, as previsões obtidas com GAM são mais precisas. 44 Anexos packages que não vinham no pacote inicial, a mboost, a doParallel e a fanplot. No entanto duas destas tinham como requisito funções de outras packages que também tiveram que ser instaladas, mas deste documento vão ser focadas apenas as três referidas. A realização das simulações e a construção de modelos em RStudio foi realizada graças à utilização da package mboost na versão 2.4-0. Este conjunto de funções permite utilizar o algoritmo do gradiente descendente para otimizar funções de custo genéricas, especificamente para este trabalho foram utilizadas as funções glmboost, gamboost e predict. Para que se possa aceder às funcionalidades da package mboost é necessária a prévia instalação de duas outras packages, nomeadamente a package stabs (Stability Selection with Error Control) e parallel que implementa funções que permitem a execução de código em paralelo em máquinas com mais do que um processador na CPU. Tal como referido no Capítulo 3 as funções glmboost e gamboost foram utilizadas para a construção dos modelos. Após a obtenção dos modelos foi necessário verificar a validade destes, para tal foi usada a função predict. Esta função recebe os dados de validação e o modelo criado pela função glmboost ou gamboost e calcula a previsão de preços com base nos argumentos que lhe são fornecidos. A grande quantidade de dados de validação para calcular, aliado a um tempo de execução superior da função predict face às funções glmboost e gamboost levaram a um tempo total para obtenção de resultados bastante alargado. Para resolver esta situação foram utilizadas as funções da package doParallel versão 1.0.8, mais especificamente as funções registerDoParallel e %dopar%. Como a máquina utilizada para este trabalho tem 8 CPU’s, mas apenas estava a utilizar um, foi possível reduzir substancialmente o tempo de cálculo das previsões com estas funções que permitiram utilizar a totalidade dos recursos disponíveis. Como requisito para a execução da doParallel foi necessário instalar as packages foreach e iterators, sendo que a primeira contém funções necessárias à criação de ciclos para computação paralela em R, e a segunda ferramentas para iterar estruturas de dados do R. Em R é possível criar gráficos usando apenas as funções das packages instaladas por defeito com o RStudio. No entanto para criar gráficos como o Gráfico 3.1 foi utilizada a package fanplot versão 3.3 que contém uma coleção de funções em R para criar gráficos de distribuições sequenciais como previsões probabilísticas ou resultados de simulações. Foi esta função que permitiu representar os quantis de previsões em bandas sombreadas de acordo com a incerteza associada aos quantis, o que permite uma melhor compreensão dos dados representados do que os gráficos de pontos e linhas que se podem criar com as funções instaladas por defeito. Package mboost No cerne de todo o trabalho esteve o package mboost e as suas funções, glmboost, gamboost e predict. De seguida será descrita e exemplificada a sua utilização. A função glmboost deve ser chamada utilizando o protótipo: 45 Software glmboost(formula, data, family, control) No qual formula é uma descrição simbólica do modelo que se pretende adaptar, data é o conjunto de dados que contém as variáveis do modelo, family é a família que já foi descrita no capítulo 3 deste documento e control é uma lista de parâmetros que controlam o algoritmo. A seguir encontra-se um exemplo retirado do script do glmboost: glm50<-glmboost(PRECO_PT~ previsao.carga+ previsao.eolica+ ADprec+ ADmod+ ADtemp+ Aprec+ Amod+ Atemp+ Lprec+ Lmod+ Ltemp+ TMprec+ TMmod+ TMtemp+ semana_anterior+ dia_anterior+ hora_antes+ duas_hora_antes, data=rel, family=QuantReg(tau=0.50), control = boost_control(mstop = 5000)) Neste exemplo é atribuída a uma variável (glm50) o modelo criado pela função glmboost, à qual foram dados como parâmetros uma descrição do modelo no qual o preço da energia é função de previsões de carga, e produção eólica, de previsões meteorológicas e de preços em instantes passados. O parâmetro family é uma regressão por quantis com τ=0.5, ou seja este é o quantil 50%. Foi também definido que o número de iterações deveria ser de 5000. A função gamboost tem um protótipo muito semelhante ao da glmboost. gamboost(formula, data, family, control) 46 Anexos Os parâmetros que são passados a esta função são do mesmo tipo dos da função glmboost. A seguir encontra-se um exemplo de chamada da função gamboost: gam95<-gamboost(PRECO_PT~ bbs(previsao.carga)+ bbs(previsao.eolica)+ bbs(ADprec)+ bbs(ADmod)+ bbs(ADtemp)+ bbs(Aprec)+ bbs(Amod)+ bbs(Atemp)+ bbs(Lprec)+ bbs(Lmod)+ bbs(Ltemp)+ bbs(TMprec)+ bbs(TMmod)+ bbs(TMtemp)+ bbs(semana_anterior)+ bbs(dia_anterior)+ bbs(hora_antes)+ bbs(duas_hora_antes)+ bbs(angulo_hora, cyclic=TRUE, boundary.knots = c(0, 2*pi))+ bbs(angulo_dds, cyclic=TRUE, boundary.knots= c(0, 2*pi))+ bbs(angulo_dia, cyclic=TRUE, boundary.knots= c(0, 2*pi))+ bbs(angulo_mes, cyclic=TRUE, boundary.knots= c(0, 2*pi)), data=rel,family=QuantReg(tau=0.95),control = boost_control(mstop = 5000)) Neste exemplo é guardado na variável gam95 o modelo criado pela função gamboost para o quantil 95%. Esta função recebe os parâmetros da mesma forma que a função glmboost, no entanto o parâmetro formula contém mais variáveis do que no exemplo do glmboost. Isto acontece porque as variáveis extra que se encontram na chamada desta função têm uma relação não linear com o preço e portanto não faria sentido incluir estas variáveis num modelo linear. A função predict pode ser usada para prever o estado da variável de resposta para novas observações, assim quando se pretendem fazer previsões com um modelo usa-se o seguinte protótipo: predict (object,newdata) 47 Software Os argumentos que a função necessita são object, que neste trabalho é o modelo criado pelo gamboost ou pelo glmboost, e newdata que são as observações das variáveis explicativas do modelo com o qual se pretende fazer previsões. A seguir é apresentado um exemplo da chamada da função: forecast50<-predict(gam50, newdata = validacao) No caso apresentado é passado à função o modelo do gamboost para o quantil 50% e os dados do conjunto de validação e o resultado da execução da função é guardado na variável forecast50. O resultado desta função é do mesmo tamanho do parâmetro newdata, ou seja como neste caso validacao é uma matriz com uma coluna por cada variável explicativa do modelo e 6022 linhas, uma por cada hora do período de validação considerado, a variável forecast50 vai ser um vetor de 6022 previsões, também uma por cada hora do período de validação. Avaliação dos Resultados Obtidos Para aferir a qualidade das previsões realizadas foram utilizadas três funções que não pertencem a nenhuma das packages consideradas mas foram fornecidas pelo Doutor Ricardo Bessa, orientador deste trabalho. Cada uma das funções produz informações quanto a três métricas sendo estas a calibração dos dados, a amplitude dos intervalos e CRPS (Continuous Rank Probability Score). A função que avalia a calibração tem o seguinte protótipo: Quantiles_Calibration(qForecast,observation) Onde qForecast é uma matriz composta pelos vetores das previsões de cada quantil e observation é a referência dos preços, ou seja os preços que efetivamente se verificaram durante o período dos dados de validação. Esta função retorna um vetor de 19 elementos, um por cada quantil. A função que avalia a amplitude dos intervalos tem o protótipo muito semelhante ao da calibração: Quantiles_Sharpness(qForecast,observation,quantiles=seq(from=5,to=95,by=5)) Apesar de no protótipo existir mais um parâmetro do que na função calibração, este já se encontra definido e portanto para a chamada destas funções basta passar os mesmos dados que 48 Anexos são passados à função calibração. Já o resultado obtido da aplicação desta função é uma matriz, com duas colunas e 9 linhas. A função CRPS tem o protótipo: Quantiles_CRPS(qForecast,observation,quantiles=seq(from=.05,to=.95,by=.05),Ymax,Ymin) Para a chamada desta função é necessário a introdução de dois novos dados, o maior valor presente na matriz das previsões e o menor valor presente nessa mesma matriz, os parâmetros Ymax e Ymin respetivamente. A execução desta função produz um vetor de 24 elementos, sendo cada um o CRPS por hora.