Séries Temporais
lsj
June 17, 2023
aula de hoje:
- O sistema binário de raios-X GX 5-1
- Curvas de luz estelares da missão Kepler
1. O sistema binário de raios-X GX 5-1
adaptação do Time series analysis, Summer School in Statistics for Astronomers, Eric D. Feigelson, 2021
https://sites.psu.edu/astrostatistics/files/2021/02/Tutorial_4_Time_Series.html
Ver também Modern Statistical Methods for Astronomy with R Applications, E.D.Feigelson & G.J. Babu (2012)
1.1 Introdução
GX 5-1 é um sistema binário que é forte emissor de raios-X. Nesse sistema, uma estrela normal transfere massa para uma estrela de neutrons via um disco de acresção, e os raios-X são emitidos na região mais interna, mais quente, desse disco.
O arquivo GX.dat contêm 65536 medidas de contagens de fotons em raios-X em bins igualmente espaçados por 1/128 segundos e foram obtidos pelo satélite japonês Ginga na década de 80.
As taxas de contagem são um processo poissoniano, mas como as contagens por bin são relativamente altas, podem ser aproximadas como um processo gaussiano.
leitura dos dados:
### Time series analysis
### Summer School in Statistics for Astronomers
### Eric D. Feigelson Spring 2021
# Read the dataset
GX.dat <- scan("GX.dat")
head(GX.dat)## [1] 79 70 68 58 50 59summary(GX.dat)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 30.00 62.00 69.00 68.82 76.00 116.00Note que os dados serão colocados num objeto de classe ‘ts’ (“time series”), adequado para a análise de séries temporais em R:
# Convert to a class 'ts' (time series) R object
# tempo: note que a série temporal corresponde a 65536/128 = 512 segundos
GX.time <- seq(from=0, to=512, length.out=length(GX.dat))
GX.ts <- ts(GX.dat, GX.time) ; GX.ts.offset <- ts(GX.dat-30, GX.time)
str(GX.ts)## Time-Series [1:65536] from -0.992 to 65534: 79 70 68 58 50 59 52 55 61 64 ...# visualização dos dados:
plot.ts(GX.ts, ylab='GX 5-1 counts', xlab='Time (x 1/128 sec)', cex.lab=1.3, cex.axis=1.3, lwd=0.5,col='gold3')
Vamos comparar a distribuição de contagens com uma distribuição normal. Preste atenção em como são calculadas a média e o desvio padrão dessa gaussiana
# opções gráficas
options(jupyter.plot_scale=1)
options(repr.plot.width = 7, repr.plot.height = 5)
hist(GX.dat, breaks=100, xlim=c(40,100), ylim=c(0,3500), xlab='GX 5-1 counts',
font=2, font.lab=2, main='')
curve(dnorm(x,mean=mean(GX.dat), sd=sqrt(mean(GX.dat)))*65536, lwd=3, add=T)
# comparação da largura dos dados com uma gaussiana
sd(GX.dat) / sqrt(mean(GX.dat)) ## [1] 1.241755Note que os dados têm uma largura 1.24 vezes maior comparada a um ruído gaussiano: aí tem coisa!
Vamos dar um zoom na série, vendo \(\sim10\)% dos dados:
plot.ts(GX.ts[1:6000], ylab='GX 5-1 counts', xlab='Time (x 1/128 sec)', cex.lab=1.3, cex.axis=1.3,col='cadetblue') 
será que tem algum comportamento semiperiódico presente?
Vamos exibir agora várias versões filtradas da série:
plot(GX.time,GX.dat, ylim=c(-10,115), xlab='Time (sec)', ylab='GX 5-1 counts',
cex.lab=1.3, cex.axis=1.3, type='n') # set up plot window but don't show any data
lines(ksmooth(GX.time, GX.dat+30, 'normal', bandwidth=7), lwd=2)
text(450, 110, 'Normal kernel') # Gaussian kernel density estimator with 7 bin FWHM bandwidth
lines(filter(GX.ts, sides=2, rep(1,7)/7), lwd=2)
text(450, 85, 'Moving average') # Moving average smoother with 7 bin bandwidth
lines(kernapply(GX.ts.offset, kernel('modified.daniell', 7)), lwd=2)
text(450, 50, 'Modified Daniell') # Moving average smoother with 1/2-weight at the end values of the span
lines(supsmu(GX.time, GX.dat-60, span=0.01), lwd=2)
text(400, 20, "Friedman's super-smoother") # A smoother with adaptive bandwidth from Friedman (1984)
lines(lowess(GX.time, GX.dat-80, 0.02), lwd=2)
text(400, 0, 'LOWESS local regression') # Cleveland's (1979) robust local polynomial smoother
Algumas das curvas acima sugerem variações com escalas de tempo de 20 a 50 segundos.
1.2 A função de autocorrelação
Uma ferramenta importante para analisar as séries temporais é a função de autocorrelação: ela analisa a correlação de uma variável entre instantes separados por um certo intervalo \(k\).
Dada uma série temporal \(\{x_1, x_2, ...\}\), a ACF(k) é a fração da variância total devido a valores com atraso de k intervalos: \[ ACF(k) = {\sum_{t=1}^{n-k} (x_t - \bar x)(x_{t-k}- \bar x) \over \sum_{t=1}^{n} (x_t - \bar x)^2} \]
A função de autocorrelação parcial (PACF) com atraso k dá a autocorrelação em k depois de remover os efeitos das correlações em atrasos menores.
par(mfrow=c(1,2))
acf(GX.dat)
pacf(GX.dat)
# as linhas tracejadas azuis representam o intervalo de confiança de 95% A autocorrelação com k=1 é maior que 20%, o que mostra uma forte dependência do valor atual de valores anteriores.
Parece haver uma componente quasiperiódica com \(P \sim 6\) segundos (k=5).
A ACF mostra que essa estrutura é significativa, excedendo o intervalo de confiança de 95% (linhas tracejadas) até \(k \sim 30\).
1.3 Periodogramas
Periodogramas são outra ferramenta importante na análise de séries temporais.
Periodograma de Schuster (clássico):
# Spectral periodogram
par(mfrow=c(1,1))
spec.pgram(GX.ts, log='no', main='')
Periodograma de Schuster suavisado:
# Spectral periodogram
par(mfrow=c(2,1))
spec.pgram(GX.ts, spans=50, log='no', main='', sub='')
spec.pgram(GX.ts, spans=200, taper=0.15, log='no', main='', sub='')
Os periodogramas mostram estruturas associadas a dois tipos de processos: o crescimento em direção a baixas frequências (abaixo de 0.05) é devido a ruído (red noise), e o pico alargado, em 0.17-0.23, é devido às oscilações quasiperiódicas.
1.4 Um modelo autoregressivo
Vamos agora considerar um modelo autoregressivo simples, onde o valor atual depende linearmente de um certo número de valores em tempos anteriores. Nesse caso, o algoritmo encontra que o melhor modelo (usando AIC) tem 27 coeficientes:
# Autoregressive modeling: AR
ARmod <- ar(GX.ts, method='ols')
print(ARmod)##
## Call:
## ar(x = GX.ts, method = "ols")
##
## Coefficients:
## 1 2 3 4 5 6 7 8
## 0.2047 0.0153 0.0046 0.0674 0.1086 0.0515 -0.0166 -0.0242
## 9 10 11 12 13 14 15 16
## 0.0059 0.0430 0.0370 -0.0053 -0.0097 -0.0073 0.0032 0.0241
## 17 18 19 20 21 22 23 24
## -0.0010 -0.0064 -0.0167 0.0040 0.0177 0.0034 0.0025 -0.0067
## 25 26 27
## -0.0050 0.0073 0.0113
##
## Intercept: -4.655e-05 (0.03832)
##
## Order selected 27 sigma^2 estimated as 96.18ARspec <- spec.ar(GX.ts, plot=F)
GXspec <- spec.pgram(GX.ts, span=101, main='', sub='', lwd=2)
lines(ARspec$freq, ARspec$spec, col='green', lwd=2)
text(0.4,450, cex=2, paste0('AR(', ARmod$order, ')'))
Uma forma de verificar se o modelo autoregressivo ajustado é bom, é calculando a função de autocorrelação dos resíduos do ajuste: espera-se que os resíduos não estejam correlacionados!
acf(na.omit(ARmod$resid), main='ACF dos resíduos')
Pode-se notar que o modelo ajusta muito bem a ACF, sem introduzir nenhuma componente periódica.
Existem vários pacotes que implementam modelos autoregressivos em R. Exemplos adicionais de aplicação a GX 5-1 podem ser encontrados em https://sites.psu.edu/astrostatistics/files/2021/02/Tutorial_4_Time_Series.html
A análise desta seção nos leva a concluir que GX 5-1 exibe variabilidade estocástica de baixa amplitude, com oscilações semiperiódicas de 20-50 segundos.
2. Curvas de luz estelares da missão Kepler
2.1 Introdução
A missão Kepler mediu com alta precisão o brilho na faixa do visível de cerca de 150 mil estrelas cada 29.4 minutos por quase 4 anos, tendo detectado milhares de trânsitos exoplanetários, onde a passagem de um planeta em frente a uma estrela provoca uma diminuição de seu brilho.
Nessa seção vamos considerar um caso de variabilidade estelar que não têm nada a ver com exoplanetas e que ilustra as dificuldades que podem aparecer na busca por esses objetos.
Esse exemplo considera séries temporais igualmente espaçadas. Mas as séries observadas têm, em geral, muitos dados faltantes e os métodos que vamos considerar exigem séries completas. A solução é imputar dados, isto é, atribuir valores aos dados faltantes (NA, Not Available) usando algum procedimento estatístico.
2.2 Curva de luz de KIC 007609553
Essa estrela mostra variabilidade quase-periódica devido à presença de rotação e manchas estelares.
# leitura dos dados: vetor pré-processado com medidas de fluxo igualmente espaçadas no tempo
Kepler2 <- read.table('Kepler2.dat')[[1]]
length(Kepler2)## [1] 71427# fração de NAs:
length(which(is.na(Kepler2))) / length(Kepler2) ## [1] 0.2883783# ~29% NAs
# visualização
par(mfrow=c(1,1))
plot(Kepler2, type='l', xlab='Time')
Note que 29% dos dados são NAs. Vamos substituir os NAs por valores “imputados”.
Há inúmeras técnicas de imputação; vamos usar um filtro de Kalman:
Em estatística, o filtro de Kalman é um método matemático criado por Rudolf Kalman. Seu propósito é utilizar medições de grandezas realizadas ao longo do tempo (contaminadas com ruído e outras incertezas) e gerar resultados que tendam a se aproximar dos valores reais das grandezas medidas e valores associados. O filtro de Kalman apresenta diversas aplicações e é uma parte essencial do desenvolvimento de tecnologias espaciais e militares. (Wikipedia)
library(imputeTS)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zooKepler2_impute <- na_kalman(Kepler2) # impute NAs with Kalman smoother## Warning in stats::StructTS(data, ...): possible convergence problem: 'optim'
## gave code = 52 and message 'ERROR: ABNORMAL_TERMINATION_IN_LNSRCH'# veja o antes e o depois da imputação:
par(mfrow=c(2,1)) ; par(mar=c(5,4,1,2))
plot(Kepler2, type='l', xlim=c(60000,70000), xlab='Time', ylab='Kepler 2')
plot(Kepler2_impute, type='l', xlim=c(60000,70000), xlab='Time',
ylab='Kepler 2 imputed')
Vamos examinar esta curva de luz (imputada) com 2 periodogramas:
- análise de Fourier usando a função spec.pgram
- o periodograma de Lomb-Scargle, muito usado em astronomia
# periodograma com análise de Fourier
par(mfrow=c(1,1))
spec.pgram(Kepler2_impute, xlim=c(0,0.005), spans=5, taper=0.0,
main='', ylab='Fourier', sub='')
# periodograma de Lomb-Scargle
library(lomb)
lsp(Kepler2, from=0.00, to=0.005, ylab='Lomb-Scargle', main='')
Ambas as técnicas sugerem dois períodos principais na curva de luz deste objeto.
Há muitas outras técnicas que todos os interessados em trabalhar com séries temporais precisam conhecer!
Exercício: análise do número de manchas solares
A série do número de manchas dolares com o tempo é bastante complexa, como pode ser vislumbrado por esse exercício.
A tabela SN_m_tot_V2.0.txt contém o número mensal de manchas solares entre 1749 e 2023 (download em 18/06/2023), como disponível em https://www.sidc.be/SILSO/infosnmtot
A tabela contém: o ano, o mês, o ano com a fração do ano no meio do mês correspondente, o número de manchas solares, o desvio padrão do número de manchas se disponível, o número de observações se disponível.
- Leia os dados da tabela e faça um gráfico do número de manchas solares em função do tempo
- Suavize essa figura usando um kernel gaussiano com largura de banda de 5 meses. Comente o resultado: como se compara com a curva do item 1?
- Faça um plot da função de autocorrelação. Como você interpreta o resultado? Atenção: para visualizar a acf considere lags de até 1000 meses (faça ?acf e use o parâmetro lag.max)
- Faça um periodograma de Lomb-Scargle. Comente sobre as frequências relevantes que você vê.
- Faça um modelo autoregressivo AR. Compare com o ajuste obtido acima para GX 5-1.