VII.8 O envelope espectral

O conceito de envelope espectral para a análise espectral e escalonamento de séries temporais categóricas foi introduzido pela primeira vez em Stoffer et al. (1993). Desde então, a ideia foi estendida em várias direções, não apenas restrita a séries temporais categóricas, e iremos explorar esses problemas também. Primeiro, damos uma breve introdução ao conceito de escalonamento de séries temporais.

O envelope espectral foi motivado por colaborações com pesquisadores que coletaram séries temporais de valores categóricos com interesse no comportamento cíclico dos dados. Por exemplo, no vetor EEG mostramos o estado de sono por minuto de uma criança obtido em um estudo sobre os efeitos da exposição pré-natal ao álcool. Detalhes podem ser encontrados em Stoffer et al. (1988), mas brevemente, um registro eletroencefalográfico (EEG) do sono de aproximadamente duas horas é obtido em um bebé 24 a 36 horas após o nascimento e o registro é avaliado por um neurologista pediátrico para o estado de sono. Existem dois tipos principais de sono: movimento ocular não rápido (NON-REM), também conhecido como sono tranquilo e movimento ocular rápido (REM), também conhecido como sono ativo. Além disso, existem quatro estágios de NON-REM (NR1-NR4), com NR1 sendo o mais ativo dos quatro estados e, finalmente, acordado (AW), que ocorre naturalmente durante a noite. Este bebê em particular nunca estava acordado durante o estudo.

Não é muito difícil notar um padrão nos dados se nos concentrarmos nos estados de sono REM versus NON-REM. Porém, seria difícil tentar avaliar padrões em uma sequência mais longa ou se houvesse mais categorias, sem algum auxílio gráfico. Um método simples seria dimensionar os dados, ou seja, atribuir valores numéricos às categorias e, em seguida, desenhar um gráfico de tempo das escalas. Como os estados têm uma ordem, uma escala óbvia seria \begin{equation} NR4 \, = \, 1, \quad NR3 \, = \, 2, \quad NR2 \, = \, 3, \quad NR1 \, = \, 4, \quad REM \, = \, 5, \quad AW \, = \, 6 \end{equation} e a Figura VII.22 mostra o gráfico de tempo usando essa escala. Outra escala interessante pode ser combinar os estados silenciosos e os estados ativos: \begin{equation} NR4 \, = \, NR3 \, = \, NR2 \, = \, NR1 \, = \, 0, \quad REM \, = \, 1, \quad AW \, = \, 2\cdot \end{equation}

O gráfico de tempo usando a segunda codificação seria semelhante à Figura VII.22 no que diz respeito ao comportamento cíclico, dentro e fora do sono tranquilo, do padrão de sono deste bebê. A Figura VII.22 mostra o periodograma dos dados do sono usando a primeira escala. Existe um grande pico na frequência correspondente a um ciclo a cada 60 minutos. Como podemos imaginar, a aparência geral do periodograma usando a segunda escala (não mostrado) é semelhante à Figura VII.22. A maioria de nós se sentiria confortável com essa análise, embora tenhamos feito uma escolha arbitrária sobre a escala específica. É evidente a partir dos dados, sem qualquer escala, que se o interesse for no ciclo do sono infantil, este estudo do sono em particular indica que o bebê faz um ciclo entre o sono ativo e o tranquilo a uma taxa de cerca de um ciclo por hora.

Como identificamos que o pico do periodograma corresponde a uma frequência de aproximadamente um ciclo a cada 60 minutos? com o seguinte código respondemos a pergunta e vemos que o pico acontece exatamente a cada 54 minutos:

A intuição usada no exemplo anterior se perde quando consideramos uma longa sequência de DNA. Resumidamente, uma fita de DNA pode ser vista como uma longa sequência de nucleotídeos ligados. Cada nucleotídeo é composto de uma base nitrogenada, um açúcar de cinco carbonos e um grupo fosfato. Existem quatro bases diferentes e podem ser agrupadas por tamanho: as pirimidinas, timina (T) e citosina (C), e as purinas, adenina (A) e guanina (G).

Os nucleotídeos são ligados entre si por uma espinha dorsal de grupos alternados de a&ccedul;úcar e fosfato com o carbono 5to. de um açúcar ligado ao carbono 3ro. do próximo, dando a direção da cadeia. As moléculas de DNA ocorrem naturalmente como uma dupla hélice composta de fitas polinucleotídicas com as bases voltadas para dentro. As duas fitas são complementares, por isso é suficiente representar uma molécula de DNA por uma sequência de bases em uma única fita. Assim, uma fita de DNA pode ser representada como uma sequência de letras, denominadas pares de bases (bp), do alfabeto finito {A, C, G, T}. A ordem dos nucleotídeos contém a informação genética específica do organismo. A expressão das informações armazenadas nessas moléculas é um processo complexo de vários estágios.

Uma tarefa importante é traduzir as informações armazenadas nas sequências codificadoras de proteínas (CDS) do DNA. Um problema comum na análise de dados de longa sequência de DNA é a identificação CDS dispersos por toda a sequência e separados por regiões de não codificação, que constituem a maior parte do DNA. Abaixo mostramos parte da sequência de DNA do vírus Epstein-Barr (EBV). Toda a sequência de DNA do EBV consiste em aproximadamente 172.000 pb.

Poderíamos tentar escalar de acordo com o alfabeto purina-pirimidina, que é A = G = 0 e C = T = 1, mas isso não é necessariamente de interesse para todos os CDS do EBV. Existem numerosos alfabetos de interesse possíveis. Por exemplo, podemos nos concentrar no alfabeto de ligações de hidrogénio forte-fraco C = G = 0 e A = T = 1. Embora os cálculos do modelo, bem como os dados experimentais concordem fortemente que algum tipo de sinal periódico existe em certas sequências de DNA, um grande desacordo sobre o tipo exato de periodicidade existe. Além disso, existe um desacordo sobre quais alfabetos de nucleotídeos estão envolvidos nos sinais.

Se considerarmos a abordagem ingênua de atribuir arbitrariamente valores numéricos (escalas) às categorias e, em seguida, prosseguir com uma análise espectral, o resultado dependerá da atribuição particular de valores numéricos. Por exemplo, considere a sequência artificial ACGTACGTACGT .... Então, definir A = G = 0 e C = T = 1 resulta na sequência numérica 010101010101..., ou um ciclo a cada dois pares de bases. Outra escala interessante seria A = 1, C = 2, G = 3 e T = 4, que resulta na sequência 123412341234... ou um ciclo a cada quatro bp. Neste exemplo, ambas as escalas, ou seja, {A, C, G, T} = {0, 1, 0, 1} e {A, C, G, T} = {1, 2, 3, 4}, dos nucleotídeos são interessantes e revelam diferentes propriedades da sequência. Portanto, não queremos nos concentrar em apenas um dimensionamento. Em vez disso, o foco deve ser localizar todas as escalas possíveis que revelem todos os recursos interessantes nos dados. Em vez de escolher valores arbitrariamente, a abordagem do envelope espectral seleciona escalas que ajudam a enfatizar qualquer recurso periódico que existe em uma série de tempo categórica de virtualmente qualquer comprimento de maneira rápida e automatizada. Além disso, a técnica pode ajudar a determinar se uma sequência é meramente uma atribuição aleatória de categorias.

O envelope espectral para séries temporais categóricas

Como uma descrição geral, o envelope espectral é uma técnica de componentes principais baseada em frequência aplicada a uma série temporal multivariada. Primeiro, vamos nos concentrar no conceito básico e seu uso na análise de séries temporais categóricas. Detalhes técnicos podem ser encontrados em Stoffer et al. (1993).

Resumidamente, ao estabelecer o envelope espectral para séries temporais categóricas, a questão básica de como descobrir componentes periódicos em séries temporais categóricas foi abordada. Isso foi realizado por meio de análise espectral não paramétrica como segue.

Seja $X_t$, $t=0,\pm 1,\pm 2,\cdots$ uma série temporal de valor categórico com espaço de estados finito $\mathcal{C}=\{c_1,c_2,\cdots,c_k\}$. Suponha que $X_t$ seja estacionária e que $p_j = P(X_t=c_j)>0$ para $j = 1,2,\cdots,k$. Para o vetor de parâmetros $\beta=(\beta_1,\beta_2,\cdots,\beta_k)^\top\in \mathbb{R}^k$, denotemos por $X_t(\beta)$ a série temporal estacionária de valor real correspondente à escala que atribui à categoria $c_j$ o valor numérico $\beta_j$, $j = 1,2,\cdots,k$. A densidade espectral de $X_t(\beta)$ será denotada por $f_{_{XX}}(\omega;\beta)$. O objetivo é encontrar $\beta$ de forma que a densidade espectral seja interessante em certo sentido e resumir a informação espectral pelo que é chamado de envelope espectral.

Em particular, $\beta$ é escolhido para maximizar a potência em cada frequência $\omega$ de interesse, em relação à potência total $\sigma^2(\beta) = \mbox{Var}\big(X_t(\beta)\big)$. Ou seja, escolhemos $\beta(\omega)$ em cada $\omega$ de interesse, de maneira que \begin{equation} \lambda(\omega) \, = \, \max_{\beta} \left( \dfrac{f_{_{XX}}(\omega;\beta)}{\sigma^2(\beta)}\right), \end{equation} sobre todo $\beta$ não proporcional a $\pmb{1}_k$, o vetor $k\times 1$ de uns.

Nota que $\lambda(\omega)$ não está definido se $\beta=a\pmb{1}_k$ para $2\in\mathbb{R}$ porque tal escala corresponde a atribuir a cada categoria o mesmo valor $a$; neste caso, $f_{_{XX}}(\omega;\beta)=0$ e $\sigma^2(\beta)=0$. O critério de otimalidade $\lambda(\omega)$ possui a propriedade desejável de ser invariante sob mudanças de localização e escala de $\beta$.

Como na maioria dos problemas de escala para dados categóricos, é útil representar as categorias em termos dos vetores unitários $u_1, u_2,\cdots,u_k$, onde $u_j$ representa o vetor $k\times 1$ com um na $j$-ésima linha e zeros em qualquer outro lugar. Definimos então uma série temporal estacionária $k$-dimensional $Y_t$ por $Y_t = u_j$ quando $X_t = c_j$.

A série temporal $X_t(\beta)$ pode ser obtida a partir da série temporal $Y_t$ pela relação $X_t(\beta)= \beta^\top Y_t$. Suponha que o processo vetorial $Y_t$ tenha uma densidade espectral contínua denotada por $f_{_{YY}}(\omega)$. Para cada $\omega$, $f_{_{YY}}(\omega)$ é uma matriz Hermitiana $k\times k$ de valor complexo.

A relação $X_t(\beta)= \beta^\top Y_t$ implica \begin{equation} f_{_{XX}}(\omega;\beta) \, = \, \beta^\top f_{_{YY}}(\omega;\beta)\beta \, = \, \beta^\top f_{_{YY}}^{^{re}}(\omega;\beta)\beta, \end{equation} onde $f_{_{YY}}^{^{re}}(\omega;\beta)$ denota a parte real de $f_{_{YY}}(\omega;\beta)$. Nesta seção, é mais conveniente escrever valores complexos na forma $z = z^{re} + iz^{im}$; que representa uma mudança da notação usada anteriormente. A parte imaginária desaparece da expressão porque é assimétrica, ou seja, $f_{_{YY}}^{^{im}}(\omega;\beta)^\top=-f_{_{YY}}^{^{im}}(\omega;\beta)$.

O critério de otimalidade pode, portanto, ser expresso como \begin{equation} \lambda(\omega) \, = \, \max_{\beta} \left( \dfrac{\beta^\top f^{^{re}}_{_{YY}}(\omega;\beta)\beta}{\beta^\top V \beta}\right), \end{equation} onde $V$ é a matriz de variância-covariância de $Y_t$. A escala resultante $\beta(\omega)$ é chamada de escala ótima.

O processo $Y_t$ é um processo multivariado de ponto e qualquer componente particular de $Y_t$ é o processo de ponto individual para o estado correspondente, por exemplo, o primeiro componente de $Y_t$ indica se o processo está no estado $c_1$ no tempo $t$. Para qualquer $t$ fixo, $Y_t$ representa uma única observação de um esquema de amostragem multinomial simples. Segue-se prontamente que $V = D-p\times p^\top$, onde $p = (p_1,\cdots,p_k)^\top$ e $D$ é a matriz $k\times k$ diagonal $D = \mbox{diag}\{p_1,\cdots,p_k\}$. Porque, por suposição, $p_j> 0$ para $j = 1,2,\cdots,k$, segue-se que $\mbox{posto}(V) = k - 1$ com o espaço nulo de $V$ sendo medido por $\pmb{1}_k$. Para qualquer matriz $k\times (k-1)$ de classificação de posto completo $Q$ cujas colunas sejam linearmente independentes de $\pmb{1}_k$, $Q^\top VQ$ será uma matriz $(k-1)\times (k-1)$ simétrica definida positiva.

Com a matriz $Q$ conforme definida anteriormente, defina $\lambda(\omega)$ como o maior autovalor da equação determinante \begin{equation} \left|Q^\top f^{^{re}}_{_{YY}}(\omega;\beta)Q \, - \, \lambda(\omega)Q^\top V Q\right| \, = \, 0 \end{equation} e seja $b(\omega)\in \mathbb{R}^{k-1}$ ser qualquer autovetor correspondente, isto é, \begin{equation} Q^\top f^{^{re}}_{_{YY}}(\omega;\beta)Q b(\omega) \, = \, \lambda(\omega)Q^\top V Q b(\omega)\cdot \end{equation} O autovalor $\lambda(\omega)>0$ não depende da escolha de $Q$. Embora o autovetor $b(\omega)$ dependa da escolha particular de $Q$, a classe de equivalência das escalas associadas com $\beta(\omega)= Qb(\omega)$ não depende de $Q$. Uma escolha conveniente de $Q$ é $Q = \big( \mbox{I}_{k-1} \, | \, 0 \big)^\top$, onde $\mbox{I}_{k-1}$ é a matriz $(k-1)\times (k-1)$ identidade e 0 é o vetor $(k-1)\times 1$ de zeros. Para esta escolha, $Q^\top f^{^{re}}_{_{YY}}(\omega;\beta)Q $ e $Q^\top VQ$ são os blocos superiores de dimensões $(k-1)\times (k-1)$ de $f^{^{re}}_{_{YY}}(\omega;\beta) $ e $V$, respectivamente. Esta escolha corresponde a definir o último componente de $\beta(\omega)$ como zero.

O próprio valor $\lambda(\omega)$ tem uma interpretação útil; especificamente, $\lambda(\omega)\mbox{d}(\omega)$ representa a maior proporção da potência total que pode ser atribuída às frequências $(\omega,\omega+\mbox{d}\omega)$ para qualquer processo em escala particular $X_t(\beta)$, com o máximo sendo alcançado pela escala $\beta(\omega)$. Devido ao seu papel central, $\lambda(\omega)$ é definido como o envelope espectral de uma série temporal categórica estacionária.

O nome envelope espectral é apropriado, pois $\lambda(\omega)$ envolve o espectro padronizado de qualquer processo escalado. Ou seja, dado qualquer $\beta$ normalizado de forma que $X_t(\beta)$ tenha potência total um, $f_{_{XX}}(\omega;\beta)\leq \lambda(\beta)$ com igualdade se, e somente se, $\beta$ for proporcional a $\beta(\omega)$.

Dadas as observações $X_t$, para $t = 1,\cdots,n$, em uma série temporal categórica, formamos o processo pontual multinomial $Y_t$, para $t = 1,\cdots,n$. Então, a teoria para estimar a densidade espectral de uma série temporal multivariada com valor real pode ser aplicada para estimar $f_{_{YY}}(\omega)$, a densidade espectral $k\times k$ de $Y_t$.

Dada uma estimativa $\widehat{f}_{_{YY}}(\omega)$ de $f_{_{YY}}(\omega)$, as estimativas $\widehat{\lambda}(\omega)$ e $\widehat{\beta}(\omega)$ do envelope espectral, $\lambda(\omega)$ e as escalas correspondentes $\beta(\omega)$, podem então ser obtidas. Detalhes sobre estimação e inferência para o envelope espectral amostral e as escalas ideais podem ser encontrados em Stoffer et al. (1993), mas o principal resultado desse artigo é o seguinte:

Se $\widehat{f}_{_{YY}}(\omega)$ é um estimador espectral consistente e se para cada $j = 1,\cdots,J$, a maior raiz de $f^{^{re}}_{_{YY}}(\omega_j)$ é distinta, então \begin{equation} \{\eta_n\big(\widehat{\lambda}(\omega_j)-\lambda(\omega_j)\big)/\lambda(\omega_j) , \, \eta_n\big(\widehat{\beta}(\omega_j)-\beta(\omega_j)\big) \, : \, j=1,\cdots,J\} \end{equation} converge, quando $n\to\infty$, conjuntamente em distribuição para distribuições normais independentes de média zero, a primeira das quais é normal padrão; a estrutura de covariância assintótica de $\widehat{\beta}(\omega)$ é discutida em Stoffer et al. (1993).

A pesquisa do pico para a estimativa de envelope espectral suavizada pode ser auxiliada usando a seguinte aproximação. Usando uma expansão de Taylor de primeira ordem, temos \begin{equation} \log\big(\widehat{\lambda}(\omega)\big) \, \approx \, \log\big(\lambda(\omega) \big)+\dfrac{\widehat{\lambda}(\omega)-\lambda(\omega)}{\lambda(\omega)}, \end{equation} então $\eta_n\big(\widehat{\lambda}(\omega_j)-\lambda(\omega_j)\big)$ é aproximadamente normal padrão. Segue-se que \begin{equation} \mbox{E}\Big( \log\big(\widehat{\lambda}(\omega)\big)\Big) \, \approx \, \log\big(\lambda(\omega)\big) \end{equation} e \begin{equation} \mbox{Var}\Big( \log\big(\widehat{\lambda}(\omega)\big)\Big) \, \approx \, \eta_n^{-2}\cdot \end{equation} Se nenhum sinal estiver presente em uma sequência de comprimento $n$, esperamos $\lambda(j/n)\approx 2/n$ para $1 < j < n/2$ e, portanto, aproximadamente $(1-\alpha)\times 100$% do tempo, $\log\big(\widehat{\lambda}(\omega)\big)$ será menor que $\log(2/n)+(z_{\alpha}/\eta_n)$, onde $z_\alpha$ é o percentil $(1-\alpha)$ da cauda superior da distribuição normal padrão. Exponenciando, o valor crítico $\alpha$ para $\widehat{\lambda}(\omega)$ torna-se $(2/n)\exp(z_{\alpha}/\eta_n)$. Os valores úteis de $z_\alpha$ são $z_{0.001} = 3.09$, $z_{0.0001} = 3.71$ e $z_{0.00001} = 4.26$ e, de acordo com nossa experiência, o limiar nesses níveis funciona bem.

Exemplo VII.16. Análise espectral de sequências de DNA.

Para ajudar a entender a metodologia, mostramos instruções explícitas para os cálculos envolvidos na estimação do envelope espectral de uma sequência de DNA $X_t$, para $t = 1,\cdots,n$, usando o alfabeto de nucleotídeos.

(i) Neste exemplo, mantemos a escala para T fixada em zero. Neste caso, formamos os vetores $3\times 1$ de dados $Y_t$: \begin{array}{cc} Y_t \, = \, (1,0,0)^\top \, \mbox{ se } X_t=\mbox{A}; & Y_t \, = \, (0,1,0)^\top \, \mbox{ se } X_t=\mbox{C}; \\ Y_t \, = \, (0,0,1)^\top \, \mbox{ se } X_t=\mbox{G}; & Y_t \, = \, (0,0,0)^\top \, \mbox{ se } X_t=\mbox{T}\cdot \end{array} O vetor de escala é $\beta=(\beta_1,\beta_2,\beta_3)^\top$ e o processo em escala é $X_t(\beta) = \beta^\top Y_t$.
(ii) Calcular a transformada discreta de Fourier (DFT) dos dados \begin{equation} Y(j/n) \, = \, \dfrac{1}{n}\sum_{t=1}^n Y_t \exp\big( -2\pi \, i \, tj/n\big)\cdot \end{equation} Observe que $Y(j/n)$ é um vetor $3\times 1$ de valor complexo. Calcule o periodograma, $I(j/n) = Y(j/n)Y^*(j/n)$ para $j=1,\cdots,[n/2]$ e guarde apenas a parte real, $I^{^{re}}(j/n)$.
(iii) Suavize o $I^{^{re}}(j/n)$ para obter uma estimativa de $f^{^{re}}(j/n)$. Sejam $\{h_k \, : \, k=0,\pm 1,\cdots, \pm m\}$ os pesos descritos na Seção IV.4. Calcular \begin{equation} \widehat{f}^{^{re}}(j/n) \, = \, \sum_{k=-m}^m h_k I^{^{re}}(j/n+k/n)\cdot \end{equation}
(iv) Calcule a matriz $3\times 3$ de variâncias-covariâncias amostral, \begin{equation} S_{_{YY}} \, = \, \dfrac{1}{n}\big(Y_t-\overline{Y}_t\big)\big(Y_t-\overline{Y}_t\big)^\top, \end{equation} onde $\overline{Y}=\frac{1}{n}\sum_{t=1}^n Y_t$ é a médias amostral dos dados.
(v) Para cada $\omega_j = j/n$, $j = 0,1,\cdots,[n/2]$, determinar o maior autovalor e o autovetor correspondente da matriz \begin{equation} \dfrac{2}{n}S_{_{YY}}^{-1/2}\widehat{f}^{^{re}}(\omega_j)S_{_{YY}}^{-1/2}\cdot \end{equation} Observe que $S_{_{YY}}^{1/2}$ é a raiz quadrada única da matriz $S_{_{YY}}$.
(vi) O envelope espectral amostral $\widehat{\lambda}(\omega_j)$ é o autovalor obtido na etapa anterior. Se $b(\omega_j)$ denota o autovetor obtido na etapa anterior, a escala amostral ótima é \begin{equation} \widehat{\beta}(\omega_j) \, = \, S_{_{YY}}^{-1/2}b(\omega_j), \end{equation} isso resultará em três valores, o valor correspondente à quarta categoria, $\mbox{T}$, sendo mantido fixo em zero.

Exemplo VII.17. Análise de um gene do vírus Epstein-Barr.

Neste exemplo, nos concentramos em uma análise dinâmica ou janela deslizante do gene BNRF1 (bp 1736-5689) do Epstein-Barr. A Figura VII.23 mostra a estimativa do envelope espectral de toda a sequência de codificação (comprimento 3954 bp). A figura também mostra um sinal forte na frequência 1/3; a escala ideal correspondente foi A = 0.10, C = 0.61, G = 0.78, T = 0, que indica que o sinal está no alfabeto de ligação forte-fraca, S = {C, G} e W = {A, T}.

A Figura VII.24 mostra o resultado do cálculo do envelope espectral sobre três janelas não sobrepostas de 1000 bp e uma janela de 954 bp, através do CDS, a saber: o primeiro, segundo, terceiro e quarto trimestres do BNRF1.

Um limite de significância aproximado de 0.0001 é 0.69%. Os primeiros três quartos contêm o sinal na frequência 1/3 (Figura VII.24 a-c); as escalas ideais amostrais correspondentes para as três primeiras janelas foram (a) A = 0.01, C = 0.71, G = 0.71, T = 0; (b) A = 0.08, C = 0.71, G = 0.70, T = 0; (c) A = 0.20, C = 0.58, G = 0.79, T = 0. As primeiras duas janelas são consistentes com a análise geral. A terceira seção, no entanto, mostra alguns pequenos desvios do alfabeto da ligação forte-fraca. O resultado mais interessante é que a quarta janela mostra que nenhum sinal est&aacte; presente. Isso leva à conjectura de que o quarto trimestre de BNRF1 do Epstein-Barr é, na verdade, não codificador.

O código R para este exemplo é o seguinte:

> library(astsa) > u = factor(bnrf1ebv) # first, input the data as factors and then > x = model.matrix(~u-1)[,1:3] # make an indicator matrix > # x = x[1:1000,] # select subsequence if desired > Var = var(x) # var-cov matrix > xspec = mvspec(x, spans=c(7,7), plot=FALSE) > fxxr = Re(xspec$fxx) # fxxr is real(fxx) > # compute Q = Var^-1/2 > ev = eigen(Var) > Q = ev$vectors%*%diag(1/sqrt(ev$values))%*%t(ev$vectors) > # compute spec envelope and scale vectors > num = xspec$n.used # sample size used for FFT > nfreq = length(xspec$freq) # number of freqs used > specenv = matrix(0,nfreq,1) # initialize the spec envelope > beta = matrix(0,nfreq,3) # initialize the scale vectors > for (k in 1:nfreq){ ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num, symmetric=TRUE) specenv[k] = ev$values[1] # spec env at freq k/n is max evalue b = Q%*%ev$vectors[,1] # beta at freq k/n beta[k,] = b/sqrt(sum(b^2)) } # helps to normalize beta > # output and graphics > frequency = xspec$freq > par(mfrow=c(1,1), mar=c(5,4,0,0), oma=c(0,0,2,2), mgp = c(2.6,0.6,0)) > plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)") > grid() > # add significance threshold to plot > m = xspec$kernel$m > # add significance threshold to plot > m = xspec$kernel$m > etainv = sqrt(sum(xspec$kernel[-m:m]^2)) > thresh=100*(2/num)*exp(qnorm(.9999)*etainv) > abline(h=thresh, lty=6, col=4) > which.max(specenv) [1] 1333 > frequency[1333] [1] 0.33325 > 1/frequency[1333] [1] 3.00075 > # details > output = cbind(frequency, specenv, beta) > colnames(output) = c("freq","specenv", "A", "C", "G") > round(output,3)[1:10,] freq specenv A C G [1,] 0.000 0.001 0.368 0.143 0.919 [2,] 0.000 0.001 0.367 0.135 0.920 [3,] 0.001 0.001 0.422 0.122 0.898 [4,] 0.001 0.001 0.575 0.136 0.807 [5,] 0.001 0.001 0.745 0.141 0.652 [6,] 0.002 0.001 0.864 0.031 0.503 [7,] 0.002 0.001 0.918 -0.199 0.343 [8,] 0.002 0.001 0.884 -0.428 0.187 [9,] 0.002 0.001 0.804 -0.592 0.056 [10,] 0.003 0.001 0.727 -0.687 -0.028

O envelope espectral para séries temporais contínuas

O conceito de envelope espectral para séries temporais categóricas foi estendido para séries temporais contínuas, $\{X_t \, : \, t = 0,\pm 1,\pm 2,\cdots\}$ em McDougall et al. (1997). O processo $X_t$ pode ter valor vetorial, mas aqui nos concentraremos no caso univariado. Mais detalhes podem ser encontrados em McDougall et al. (1997). O conceito é semelhante à busca por projeção (Friedman and Stuetzle, 1981).

Seja $\mathcal{G}$ um espaço vetorial $k$-dimensional de transformações reais contínuas com $\{g_1,\cdots,g_k\}$ sendo um conjunto de funções básicas que satisfazem $\mbox{E}\big( g_i^2(X_t)\big)<\infty$, $i = 1,\cdots,k$. Análogo ao caso da série temporais categóricas, defina a série temporal escalonada em relação ao conjunto $\mathcal{G}$ para ser o processo com valor real \begin{equation} X_t(\beta) \, = \, \beta^\top Y_t \, = \, \beta_1g_1(X_t)+\cdots \beta_k g_k(X_t) \end{equation} obtida do processo vetorial \begin{equation} Y_t \, = \, \big(g_1(X_t),\cdots,g_k(X_t) \big)^\top, \end{equation} onde $\beta=(\beta_1,\cdots,\beta_k)^\top\in\mathbb{R}^k$.

Se o processo vetorial $Y_t$, for assumido como tendo uma densidade espectral contínua $f_{_{YY}}(\omega)$, então $X_t(\beta)$ terá uma densidade espectral contínua $f_{_{XX}}(\omega;\beta)$ para todo $\beta\neq 0$. Observando que \begin{equation} f_{_{XX}}(\omega;\beta) \, = \, \beta^\top f_{_{YY}}(\omega)\beta \, = \, \beta^\top f^{^{re}}(\omega)\beta \end{equation} e \begin{equation} \sigma^2(\beta) \, = \, \mbox{Var}\big( X_t(\beta)\big) \, = \, \beta^\top V \beta, \end{equation} onde $V=\mbox{Var}(Y_t)$ é assumida como definida positivo; o critério de otimalidade \begin{equation} \lambda(\omega) \, = \, \sup_{\beta\neq 0} \left( \dfrac{\beta^\top f^{^{re}}(\omega)\beta}{\beta^\top V \beta}\right), \end{equation} está bem definido e representa a maior proporção da potência total que pode ser atribuída à frequência $\omega$ para qualquer processo em escala particular $X_t(\beta)$. Esta interpretação de $\lambda(\omega)$ é consistente com a noção do envelope espectral introduzido na seção anterior e fornece a seguinte definição de trabalho: O envelope espectral de uma série temporal com relação ao espaço $\mathcal{G}$ é definido como $\lambda(\omega)$.

A solução para este problema, como no caso categórico, é alcançada encontrando-se o maior escalar $\lambda(\omega)$ tal que \begin{equation} f^{^{re}}(\omega)\beta(\omega) \, = \, \lambda(\omega) V \beta(\omega), \end{equation} para $\beta(\omega)\neq 0$. Ou seja, $\lambda(\omega)$ é o maior autovalor de $f^{^{re}}(\omega)$ na métrica de $V$ e a escala ótima $\beta(\omega)$, é o autovetor correspondente.

Se $X_t$ é uma série temporal categórica assumindo valores no espaço de estados finito $\mathcal{S} = \{c_1,c_2,\cdots,c_k\}$, onde $c_j$ representa uma categoria particular, uma escolha apropriada para $\mathcal{G}$ é o conjunto de funções indicadoras $g_j(X_t) = I(X_t = c_j)$. Portanto, esta é uma generalização natural do caso categórico.

No caso categórico, $\mathcal{G}$ não consiste em funções $g$ linearmente independentes, mas foi fácil superar esse problema reduzindo a dimensão em um. No caso de valor vetorial, $X_t = (X_{1t},\cdots,X_{pt})^\top$, consideramos $\mathcal{G}$ como a classe de transformações de $\mathbb{R}^p$ em $\mathbb{R}$ tal que existe a densidade espectral de $g(X_t)$.

Uma classe de transformações de interesse são combinações lineares de $X_t$. Em Tiao et al. (1993), por exemplo, transformações lineares desse tipo são usadas em uma abordagem no domínio do tempo para investigar relações contemporâneas entre os componentes de séries temporais multivariadas. A estimação e a inferência para o caso real são análogas aos métodos descritos na seção anterior para o caso categórico. Consideramos um exemplo aqui; numerosos outros exemplos podem ser encontrados em McDougall et al. (1997).

Exemplo VII.18. Transformações ótimas para dados financeiros: retornos da NYSE.

Em muitas aplicações financeiras, normalmente se trata da análise dos retornos quadrados, como foi feito na Seção V.3 e na Seção VI.11. No entanto, pode haver outras transformações que fornecem mais informações do que simplesmente quadrar os dados. Por exemplo, Ding et al. (1993) que aplicou transformações da forma $|X_t|^d$, para $d\in (0,3]$, à série do mercado de ações S&P 500. Eles descobriram que a transformação de poder do retorno absoluto tem uma autocorrelação bastante alta para longas defasagens e essa propriedade é mais forte quando $d$ está em torno de 1. Eles concluíram que o "resultado parece argumentar contra as especificações do tipo ARCH com base em retornos quadrados."

Neste exemplo, examinamos os retornos da NYSE (nyse). Usamos com o conjunto gerador $\mathcal{G}=\{x, |x|,x^2\}$, que parece natural para esta análise para estimar o envelope espectral e o resultado é mostrado na Figura VII.25. Embora os dados sejam um ruído branco, eles claramente não são independentes identicamente distribuídos e uma potência considerável está presente nas frequências baixas. A presença de potência espectral em frequências muito baixas em séries econômicas destendidas tem sido freqüentemente relatada e está tipicamente associada à dependência de longo alcance.

A transformação ótima estimada perto da frequência zero $\omega = 0.001$; foi $\widehat{\beta}(0.001) = (-1,921,-2596)^\top$, o que leva à transformação \begin{equation} g(x) \, = \, -x +921 |x|-2596 x^2\cdot \end{equation} Essa transformação é mostrada na Figura VII.26. A transformação dada acima é basicamente o valor absoluto, com alguma ligeira curvatura e assimetria, para a maioria dos valores, mas o efeito dos extremos é amortecido.

O seguinte código R foi usado neste exemplo.

> library(astsa) > u = astsa::nyse # accept no substitutes > x = cbind(u, abs(u), u^2) > Var = var(x) # var-cov matrix > xspec = mvspec(x, spans=c(5,3), taper=.5, plot=FALSE) > fxxr = Re(xspec$fxx) # fxxr is real(fxx) > # compute Q = Var^-1/2 > ev = eigen(Var) > Q = ev$vectors%*%diag(1/sqrt(ev$values))%*%t(ev$vectors) > # compute spec env and scale vectors > num = xspec$n.used # sample size used for FFT > nfreq = length(xspec$freq) # number of freqs used > specenv = matrix(0,nfreq,1) # initialize the spec envelope > beta = matrix(0,nfreq,3) # initialize the scale vectors > for (k in 1:nfreq){ + ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num) # get evalues of normalized spectral matrix at freq k/n + specenv[k] = ev$values[1] # spec env at freq k/n is max evalue + b = Q%*%ev$vectors[,1] # beta at freq k/n + beta[k,] = b/b[1] # first coef is always 1 + } > # output and graphics > par(mar=c(2.5,2.75,.5,.5), mgp=c(1.5,.6,0)) > frequency = xspec$freq > plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)") > grid() > m = xspec$kernel$m > etainv = sqrt(sum(xspec$kernel[-m:m]^2)) > thresh = 100*(2/num)*exp(qnorm(.9999)*etainv)*matrix(1,nfreq,1) > lines(frequency, thresh, lty=2, col=4) > # details > b = sign(b[2])*output[2,3:5] # sign of |x| positive for beauty > output = cbind(frequency, specenv, beta) > colnames(output)=c("freq","specenv","x", "|x|", "x^2"); round(output, 4) freq specenv x |x| x^2 [1,] 0.0005 0.0109 1 58.3066 -163.2372 [2,] 0.0010 0.0115 1 -920.7998 2596.1907 [3,] 0.0015 0.0114 1 -80.4564 216.9281 [4,] 0.0020 0.0099 1 -57.9771 135.0177 [5,] 0.0025 0.0077 1 -59.4880 94.5836 [6,] 0.0030 0.0061 1 -80.1222 50.6889 [7,] 0.0035 0.0053 1 -207.6324 0.5130 [8,] 0.0040 0.0051 1 103.3062 -0.3826 [9,] 0.0045 0.0054 1 85.9370 -54.7112 [10,] 0.0050 0.0059 1 -123.8367 162.5133 [11,] 0.0055 0.0062 1 -24.0056 41.6873 [12,] 0.0060 0.0063 1 -12.1220 23.6645 [13,] 0.0065 0.0059 1 -8.4317 15.5591 [14,] 0.0070 0.0051 1 -6.5134 8.3303 [15,] 0.0075 0.0043 1 -4.9093 1.2865 [16,] 0.0080 0.0035 1 -3.6974 -5.3222 [17,] 0.0085 0.0030 1 -2.9950 -12.8396 [18,] 0.0090 0.0025 1 -2.8198 -23.6750 [19,] 0.0095 0.0022 1 -3.0274 -48.2784 [20,] 0.0100 0.0020 1 -3.4472 -125.9103 [21,] 0.0105 0.0018 1 -8.8588 -1934.0887 [22,] 0.0110 0.0018 1 -2.2095 162.6570 [23,] 0.0115 0.0018 1 -2.2436 72.1871 [24,] 0.0120 0.0020 1 -2.1775 49.3826 [25,] 0.0125 0.0022 1 -2.0621 43.3677 [26,] 0.0130 0.0021 1 -1.8706 45.0066 [27,] 0.0135 0.0020 1 -1.2180 67.7015 [28,] 0.0140 0.0021 1 -10.5157 -357.4853 [29,] 0.0145 0.0024 1 -3.7698 -48.6924 [30,] 0.0150 0.0026 1 -3.9351 -38.9165 [31,] 0.0155 0.0026 1 -6.1528 -73.6160 [32,] 0.0160 0.0025 1 18.4480 337.1603 [33,] 0.0165 0.0025 1 3.4514 82.3713 [34,] 0.0170 0.0026 1 3.0569 70.9472 [35,] 0.0175 0.0027 1 3.8469 79.5734 [36,] 0.0180 0.0026 1 4.0414 86.0489 [37,] 0.0185 0.0026 1 2.4086 63.8614 [38,] 0.0190 0.0025 1 1.1670 40.4263 [39,] 0.0195 0.0026 1 0.8364 30.1248 [40,] 0.0200 0.0027 1 1.0687 26.9573 [41,] 0.0205 0.0027 1 1.7445 28.5891 [42,] 0.0210 0.0026 1 2.9385 38.9515 [43,] 0.0215 0.0025 1 5.1067 71.4802 [44,] 0.0220 0.0023 1 6.0262 112.4082 [45,] 0.0225 0.0021 1 3.0352 92.5015 [46,] 0.0230 0.0020 1 1.3219 68.7585 [47,] 0.0235 0.0020 1 0.9257 53.4277 [48,] 0.0240 0.0022 1 1.1417 42.6927 [49,] 0.0245 0.0023 1 1.4342 35.3365 [50,] 0.0250 0.0024 1 1.4736 30.5584 [51,] 0.0255 0.0023 1 1.4025 29.1581 [52,] 0.0260 0.0022 1 1.4188 32.5456 [53,] 0.0265 0.0019 1 2.0094 45.3885 [54,] 0.0270 0.0019 1 13.9894 180.9517 [55,] 0.0275 0.0021 1 -4.9480 -27.5422 [56,] 0.0280 0.0024 1 -2.5678 -8.4593 [57,] 0.0285 0.0026 1 -1.6090 -5.3089 [58,] 0.0290 0.0026 1 -1.0982 -5.5530 [59,] 0.0295 0.0024 1 -0.9229 -8.5885 [60,] 0.0300 0.0021 1 -1.2360 -18.5742 [61,] 0.0305 0.0020 1 -3.3572 -63.8594 [62,] 0.0310 0.0021 1 12.0278 263.9838 [63,] 0.0315 0.0022 1 2.7066 70.0831 [64,] 0.0320 0.0022 1 1.8453 54.8296 [65,] 0.0325 0.0021 1 2.7576 83.7008 [66,] 0.0330 0.0021 1 -4.5705 -136.6243 [67,] 0.0335 0.0021 1 -1.8349 -52.3344 [68,] 0.0340 0.0021 1 -2.3309 -51.6034 [69,] 0.0345 0.0022 1 -3.6010 -58.1728 [70,] 0.0350 0.0023 1 -5.2708 -67.7980 [71,] 0.0355 0.0023 1 -6.3367 -75.3293 [72,] 0.0360 0.0023 1 -5.3295 -69.3717 [73,] 0.0365 0.0022 1 -3.1051 -50.6968 [74,] 0.0370 0.0021 1 -1.2742 -32.7967 [75,] 0.0375 0.0021 1 -0.3844 -23.0637 [76,] 0.0380 0.0021 1 -0.2911 -20.2568 [77,] 0.0385 0.0021 1 -0.8253 -26.4093 [78,] 0.0390 0.0020 1 -4.7325 -102.9859 [79,] 0.0395 0.0020 1 3.0850 77.0712 [80,] 0.0400 0.0020 1 0.7622 41.4081 [81,] 0.0405 0.0019 1 -0.3177 37.6300 [82,] 0.0410 0.0018 1 -1.2768 44.5158 [83,] 0.0415 0.0017 1 -2.6761 67.1266 [84,] 0.0420 0.0016 1 -6.1982 142.3070 [85,] 0.0425 0.0015 1 -30.6465 708.1127 [86,] 0.0430 0.0015 1 -81.2198 1877.2226 [87,] 0.0435 0.0015 1 -7.0923 164.7465 [88,] 0.0440 0.0014 1 -2.5703 64.0244 [89,] 0.0445 0.0014 1 -1.5715 44.8165 [90,] 0.0450 0.0014 1 -1.2862 44.8834 [91,] 0.0455 0.0013 1 -1.2262 58.2956 [92,] 0.0460 0.0013 1 -1.5795 133.7781 [93,] 0.0465 0.0013 1 0.0915 -92.7764 [94,] 0.0470 0.0013 1 -0.1651 -30.3083 [95,] 0.0475 0.0013 1 0.0137 -25.0783 [96,] 0.0480 0.0013 1 0.5280 -36.7888 [97,] 0.0485 0.0013 1 2.7786 -108.3667 [98,] 0.0490 0.0013 1 -5.8857 196.5450 [99,] 0.0495 0.0014 1 -1.9722 66.9382 [100,] 0.0500 0.0015 1 -1.2945 46.5127 [101,] 0.0505 0.0016 1 -1.0139 39.3986 [102,] 0.0510 0.0016 1 -0.8462 37.5058 [103,] 0.0515 0.0015 1 -0.6748 38.3142 [104,] 0.0520 0.0015 1 -0.4599 40.8547 [105,] 0.0525 0.0014 1 -0.3298 47.1095 [106,] 0.0530 0.0014 1 -0.5955 64.9047 [107,] 0.0535 0.0014 1 -2.7231 146.1488 [108,] 0.0540 0.0014 1 22.2657 -803.4113 [109,] 0.0545 0.0014 1 3.5619 -148.8929 [110,] 0.0550 0.0014 1 -0.8347 -79.5021 [111,] 0.0555 0.0016 1 -5.2367 -62.8310 [112,] 0.0560 0.0017 1 -29.0721 -252.3916 [113,] 0.0565 0.0017 1 5.7685 72.1579 [114,] 0.0570 0.0016 1 1.5023 42.1321 [115,] 0.0575 0.0016 1 0.5875 37.8377 [116,] 0.0580 0.0016 1 0.5732 35.9770 [117,] 0.0585 0.0016 1 0.8303 32.1439 [118,] 0.0590 0.0016 1 0.8361 28.1165 [119,] 0.0595 0.0015 1 0.5122 26.5466 [120,] 0.0600 0.0014 1 0.0715 27.2215 [121,] 0.0605 0.0014 1 -0.2978 29.1838 [122,] 0.0610 0.0012 1 -0.5597 32.5658 [123,] 0.0615 0.0011 1 -0.8296 42.1841 [124,] 0.0620 0.0010 1 -5.3360 334.2080 [125,] 0.0625 0.0011 1 -0.2360 -18.3163 [126,] 0.0630 0.0015 1 -0.5755 -5.0040 [127,] 0.0635 0.0020 1 -0.6595 -1.8800 [128,] 0.0640 0.0023 1 -0.6553 -1.0508 [129,] 0.0645 0.0022 1 -0.5503 -1.6518 [130,] 0.0650 0.0017 1 -0.2608 -5.3179 [131,] 0.0655 0.0013 1 0.9891 -25.7525 [132,] 0.0660 0.0013 1 5.1276 -92.9984 [133,] 0.0665 0.0013 1 5.4701 -94.3737 [134,] 0.0670 0.0013 1 3.3041 -58.6000 [135,] 0.0675 0.0013 1 2.3671 -45.7580 [136,] 0.0680 0.0012 1 2.0993 -47.3031 [137,] 0.0685 0.0011 1 2.7259 -79.1047 [138,] 0.0690 0.0011 1 -9.8975 387.3509 [139,] 0.0695 0.0012 1 -1.1604 59.7902 [140,] 0.0700 0.0012 1 -0.6300 40.2629 [141,] 0.0705 0.0013 1 -0.7980 41.6549 [142,] 0.0710 0.0013 1 -3.2936 94.9797 [143,] 0.0715 0.0014 1 5.1834 -95.0298 [144,] 0.0720 0.0015 1 2.8018 -45.8588 [145,] 0.0725 0.0015 1 1.8412 -33.2568 [146,] 0.0730 0.0015 1 0.4055 -18.1649 [147,] 0.0735 0.0015 1 -0.6083 -7.3157 [148,] 0.0740 0.0015 1 -0.9557 -5.4420 [149,] 0.0745 0.0013 1 -1.2919 -13.0049 [150,] 0.0750 0.0012 1 2.0491 77.1723 [151,] 0.0755 0.0014 1 0.1028 25.3431 [152,] 0.0760 0.0015 1 -0.1198 21.9134 [153,] 0.0765 0.0015 1 -0.1620 21.0437 [154,] 0.0770 0.0014 1 0.0374 20.0136 [155,] 0.0775 0.0014 1 0.4293 18.3451 [156,] 0.0780 0.0014 1 0.7117 17.2019 [157,] 0.0785 0.0013 1 0.6856 17.9673 [158,] 0.0790 0.0010 1 -0.2441 28.0712 [159,] 0.0795 0.0011 1 -4.1748 68.9244 [160,] 0.0800 0.0012 1 -5.7772 83.8180 [161,] 0.0805 0.0011 1 -5.9921 89.7256 [162,] 0.0810 0.0009 1 -10.0665 178.8503 [163,] 0.0815 0.0009 1 0.9836 -31.4059 [164,] 0.0820 0.0010 1 0.2216 -13.6028 [165,] 0.0825 0.0010 1 0.4214 -15.8080 [166,] 0.0830 0.0009 1 1.4671 -34.1289 [167,] 0.0835 0.0009 1 10.4728 -204.4560 [168,] 0.0840 0.0008 1 -6.4715 127.2476 [169,] 0.0845 0.0009 1 -2.3902 51.2155 [170,] 0.0850 0.0010 1 -1.3721 35.5850 [171,] 0.0855 0.0011 1 -0.9422 30.4892 [172,] 0.0860 0.0011 1 -0.9590 27.7353 [173,] 0.0865 0.0013 1 -1.3903 26.1408 [174,] 0.0870 0.0014 1 -1.7823 26.5196 [175,] 0.0875 0.0014 1 -2.1061 29.7715 [176,] 0.0880 0.0012 1 -2.9182 43.7514 [177,] 0.0885 0.0011 1 7.0673 -126.8521 [178,] 0.0890 0.0011 1 1.1258 -27.0295 [179,] 0.0895 0.0012 1 0.4054 -19.3034 [180,] 0.0900 0.0012 1 0.0773 -19.8375 [181,] 0.0905 0.0011 1 -0.0006 -23.4900 [182,] 0.0910 0.0011 1 0.1504 -23.9949 [183,] 0.0915 0.0011 1 0.2644 -22.4597 [184,] 0.0920 0.0011 1 0.0638 -20.9592 [185,] 0.0925 0.0011 1 -0.3460 -22.9910 [186,] 0.0930 0.0010 1 -0.7826 -44.8794 [187,] 0.0935 0.0009 1 -0.0607 266.1194 [188,] 0.0940 0.0009 1 -0.8446 47.4377 [189,] 0.0945 0.0010 1 -0.8162 33.3766 [190,] 0.0950 0.0011 1 -0.7146 29.5971 [191,] 0.0955 0.0010 1 -0.6966 29.9313 [192,] 0.0960 0.0010 1 -1.0579 37.5716 [193,] 0.0965 0.0009 1 -3.4709 84.4991 [194,] 0.0970 0.0009 1 -42.5099 819.0196 [195,] 0.0975 0.0010 1 -12.8065 235.1265 [196,] 0.0980 0.0010 1 -2.5949 53.8946 [197,] 0.0985 0.0010 1 -0.9678 28.0556 [198,] 0.0990 0.0011 1 -0.4783 22.1974 [199,] 0.0995 0.0011 1 -0.2524 20.2291 [200,] 0.1000 0.0013 1 -0.2038 18.4697 [ reached getOption("max.print") -- omitted 800 rows ] > par(mar=c(2.5,2.5,.5,.5), mgp=c(1.5,.6,0)) > # plot transform > g = function(x) { b[1]*x+b[2]*abs(x)+b[3]*x^2 } > curve(g, -.2, .2, panel.first=grid(lty=2)) > g2 = function(x) { b[2]*abs(x) } # corresponding |x| > curve(g2, -.2,.2, add=TRUE, lty=6, col=4)