Um exemplo real

Neste exemplo, vamos buscar uma curva que descreva bem o comportamento do nível dos reservatórios do sistema Cantareira, um dos principais reservatórios de água que abastecem a região metropolitana de São Paulo.

Os dados foram obtidos a partir do Portal dos Mananciais da Sabesp, ao longo de vários anos. Neste portal, as informações são atualizadas diariamente. Em minha página publico um acompanhamento desses dados, atualizado regularmente. Abaixo, reproduzo a imagem gerada no dia 01/12/2020.

Usando a informação desse portal, foi criado um arquivo de texto contendo leituras desde 2004 até o dezembro de 2020. Cada leitura está em uma linha e contém quatro colunas numéricas:

1. Dia da leitura (valor inteiro entre 1 e 31)
2. Mês da leitura (valor inteiro entre 1 e 12)
3. Ano da leitura (valor inteiro entre 2004 e 2020)
4. Percentual ocupado do volume do reservatório (valor real positivo)

Faça o download do arquivo de dados para poder acompanhar o desenvolvimento e fazer seus próprios experimentos.

Carregando e entendendo os dados

A primeira tarefa é inspecionar o conjunto de dados e conseguir exibi-lo. Para isso, vamos iniciar carregando os dados com o Octave.

load "cantareira.dat"
whos
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        C        2708x4                      86656  double

Total is 10832 elements using 86656 bytes

Com isso já descobrimos que o dado é composto de 2708 leituras. Logo exibi-las em uma tabela é impraticável e pouco informativo. Seria muito mais útil criar um gráfico.

plot(C(:,4))

Esse gráfico foi construído sem explicitar qual vetor representa as ordenadas. Neste caso, o Octave utiliza com ordenada o índice da amostra, 1 a 2708. Já podemos tirar algumas conclusões comparando este gráfico com o exibido no topo da página. Claramente, as primeira 500 amostras devem ser anteriores a 2015. As amostras seguintes (aproximadamente 2200) correspondem a um período menor (5 anos). Portanto, a amostragem do dado não é regular, ou seja as amostras contidas no dado não foram coletadas a intervalos regulares de tempo. No início da série histórica temos poucas amostras, mais espaçadas no tempo, e a partir de 2015 a amostragem deve ficou mais densa. Para exibir o gráfico corretamente precisamos utilizar como vetor de ordenadas a informação correta da data de registro de cada amostra.

Para ajustarmos uma função aos dados precisamos representar a variável dependente (porcentagem do volume do reservatório, o vetor $y$ do problema de ajuste) em função de uma variável independente e contínua. Logo, a data (dia/mês/ano) terá quer ser mapeada em uma quantidade real (será o vetor $x,$ no problema de ajuste). Não há uma única maneira de fazer isso. Podemos contar os dias corridos a partir da primeira leitura, por exemplo. Observe as primeira 4 entradas da matriz de dados.

C(1:4,:)
ans =

      1.0000      9.0000   2004.0000     51.6504
     15.0000      9.0000   2004.0000     49.5617
      1.0000     10.0000   2004.0000     47.2409
     15.0000     10.0000   2004.0000     46.0805

Vemos que a primeira leitura foi feita em 01/09/2004. Esse seria o dia 1. A leitura seguinte foi no feita em 15/09/2004, o que corresponderia ao dia 15. Já a terceira e quarta leituras corresponderiam aos dias corridos 31 e 46, respectivamente. Observe também que

C(end,:)
ans =

      1.0000     12.0000   2020.0000     47.3183

A última leitura foi tomada em 01/12/2020, que em termos de dias corridos corresponde ao dia 5935. Converter datas em dias corridos a partir de uma data inicial é algo chato de fazer, por conta da diferença de dias em meses e anos bissextos. Além disso, trabalhar com esse intervalo grande para $x$ em geral é ruim do ponto de vista numérico. Lembre que vamos precisar avaliar as funções $\phi_j$ do espaço de aproximação em todos os valores de $x.$ Isso pode levar a variações muito grandes de ordem de magnitude nas entradas da matriz $\mathsf{\Phi}.$

Ao invés de pensar em dias corridos, podemos indexar as leituras em termos do ano contínuo, ou seja, uma amostra tomada em 01/01/2010 corresponderia a $x=2010$ e uma leitura tomada em 01/02/2010 corresponderia a $x=2010 + 31/365.25 = 2010.08$ e uma amostra tomada em 01/07/2010 corresponderia a $x=2010+181/365.25 = 2010.50.$ Claro que ainda precisamos saber converter datas em dias corridos para poder computar a razão dias corridos/$365.25.$ A graça é que podemos tomar um valor médio para o número de dias em um mês ($365.25/12$) pois o pequeno erro nessa conta não chega a ser expressivo para o nosso propósito. Ainda assim, teríamos valores grandes em ordem de magnitude para a variável $x.$ Ao invés de considerar apenas o ano corrido podemos pensar no ano corrido em relação a um ano de referência. No caso do conjunto de dados, 2012 (média dos anos observados) é uma boa escolha. Com isso, ficamos com a seguinte transformação para mapear $D/M/A$ em $x,$ \begin{equation} x = \left[A + {(M-1) (365.25/12) + (D-1) \over 365.25} \right] - 2012, \end{equation} onde $365.25$ é o número médio de dias em um ano (para dar conta de anos bissextos) e $(365.25/12)$ é o número de médio de dias em um mês. Executando isto no Octave, temos:

x = ( C(:,3) + ((C(:,2)-1)*(365.25/12)+C(:,1))/365.25 ) - 2012;
[min(x) max(x)]
ans =

  -7.3306   8.9194

y = C(:,4);
plot(x+2012,y)

Agora o vetor $x$ varia em um intervalo razoável $[-7.3306,8.9194].$ Repare que no momento de gerar o gráfico adicionamos novamente 2012 ao vetor $x$ apenas para que pudéssemos ver no eixo horizontal o valor real do ano.

Apesar de totalmente negligenciada nos textos clássicos sobre Cálculo Numérico, quando trabalhamos com dados reais, esta tarefa de carregar, inspecionar e preparar os dados é fundamental para a obtenção de um bom resultado.

Resolução do problema

Com os dados preparados, passamos ao problema de ajuste, em si. O primeiro a fazer é a escolha do espaço de aproximação, ou seja, das funções $\phi_j$ que combinaremos para gerar a função de ajuste $\phi.$

A função $\phi_1(x) = 1$ deve sempre ser utilizada. Ela permitirá que o gráfico de $\phi$ desloque-se verticalmente. Observando o último gráfico acima, claramente vemos que há uma oscilação sazonal do nível da água com o período de um ano. Isso nos motiva a usar funções seno e cosseno. Porém não basta acrescentar as funções $\sin(x)$ e $\cos(x).$ Essas duas funções têm período $2\pi$. Se o período alvo é $T,$ precisamos de funções do tipo $$ \sin\left({2\pi\over T}x\right), \quad \cos\left({2\pi\over T}x\right). $$ Para oscilações com período anual, $T=1,$ poderíamos tentar usar $\sin(2\pi x)$ e $\cos(2\pi x).$

Usando como base para o espaço de aproximação $\{1,\sin(2\pi x),\cos(2\pi x)\},$ precisamos construir a matriz $\mathsf{\Phi}$ e resolver o sistemas sobredeterminado $\mathsf{\Phi}c=y.$ Até o momento só vimos como resolver esse sistema, no sentido de quadrados mínimos, através da construção do sistema normal. Em uma aula mais adiante, veremos outra forma de resolvê-lo. Por ora, saiba que o Octave consegue resolvê-lo diretamente, no sentido de quadrados mínimos, que é o que usaremos a seguir. Vejamos como fazer isso no Octave.

m = 2708;
P = [ones(m,1) sin(2*pi*x) cos(2*pi*x)];
c = P\y
c =

   50.3609
    5.3851
   -3.9387

phi = P*c;
plot(x+2012,y,'k', x+2012,phi,'r')

Claramente, conseguimos captar tanto o valor médio (por volta de 50), quanto as oscilações anuais, mas ainda estamos longe de ajustar bem a função. Certamente precisaremos de mais graus de liberdade para ajustar o comportamento complexo observado. É difícil antever quais exatamente serão as funções de base que vão dar origem a um bom ajuste. Neste ponto vamos nos permitir alguma liberdade em acrescentar mais termos, considerando bases $$ \{1,\sin(2\pi x),\cos(2\pi x), \ldots, \sin(2N\pi x),\cos(2N\pi x)\}, $$ para algum $N$ a ser testado.

N = 5;
P = ones(m,2*N+1);
for k=1:N
P(:,[2*k 2*k+1]) = [sin(2*k*pi*x) cos(2*k*pi*x)];
endfor
c = P\y;
phi = P*c;
plot(x,y,x,phi)

Ao aumentar o quantidade de termos do tipo $\sin(2k\pi x)$ e $\cos(2k\pi x)$ estamos conseguindo ajustar melhor o comportamento da curva dentro do período de um ano, mas ainda não conseguimos alterar o comportamento global. Isto porque essas funções correspondem a períodos de oscilação cada vez menores, ou períodos maiores. Aproveitando a matriz $\mathsf{\Phi}$ ($P$ no Octave), que já temos, vamos apenas acrescentar novas colunas referentes às novas funções de base.

whos P
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        P        2708x11                    238304  double

Total is 29788 elements using 238304 bytes

P(:,12:13) = [sin(x) cos(x)];
c = P\y;
phi = P*c;
plot(x,y,x,phi)

Melhorou. Vamos adicionar alguns monômios, na esperança de que isto ajude a representar esse comportamento global. Neste último teste consideramos $$ \{1,x,x^2,x^3,\sin x,\cos x, \sin(2\pi x),\cos(2\pi x),\ldots,\sin(2\cdot5\pi x),\cos(2\cdot 5\pi x)\}$$ como a seguinte base para o espaço de aproximação. O resultado está na figura a seguir.

Neste exemplo, fizemos considerações sobre período de oscilação observado no dado e usamos isso para guiar a escolha das funções de base. Talvez fosse mais adequado utilizar outra ferramenta como a transformada de Fourier discreta, que naturalmente lida com a decomposição do dado em funções trigonométricas de diferentes períodos. Mas isto não seria imediato neste problema, por conta da amostragem não regular. Além disso, uma vantagem da técnica de ajuste é que podemos mesclar funções distintas. O resultado da última figura foi obtido combinando-se 16 funções, entre monômios e funções trigonométricas. Será que apenas 16 termos da transformada de Fourier discreta produziriam resultado semelhante?

Em resumo, o problema de ajuste de curvas reside em dois aspectos fundamentais, compreender e preparar seu dado, e fazer boas escolhas para a base do espaço de aproximação.

1. Faça o download do dado usado nessa aula e tente reproduzir o que foi feito. Depois, experimente outras bases para o espaço de aproximação. Tente ajustar uma curva que represente apenas o comportamento médio dos dados, sem se preocupar em descrever oscilações anuais.

2. Tente determinar uma função que ajuste bem o comportamento da cotação do dólar. Faça o download do arquivo de dados. Neste arquivo há o registro da cotação do dólar desde 01/02/1999, quando entrou em vigor o sistema de câmbio flutuante no Brasil, até meados de outubro de 2021, momento em que o país amargava uma situação lamentável.