Análise de Série Temporal Aplicada à Medicina
em R
08 agosto 2025 20:40h
Bastão de Asclépio & Distribuição Normal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))options(warn=-1)
suppressMessages(library(astsa, warn.conflicts=FALSE))
suppressMessages(library(aTSA, warn.conflicts=FALSE))
suppressMessages(library(broom, warn.conflicts=FALSE))
suppressMessages(library(changepoint, warn.conflicts=FALSE))
suppressMessages(library(DescTools, warn.conflicts=FALSE))
suppressMessages(library(dplyr, warn.conflicts=FALSE))
suppressMessages(library(estimatr, warn.conflicts=FALSE))
suppressMessages(library(forecast, warn.conflicts=FALSE))
suppressMessages(library(ggfortify, warn.conflicts=FALSE))
suppressMessages(library(ggplot2, warn.conflicts=FALSE))
suppressMessages(library(KernSmooth, warn.conflicts=FALSE))
suppressMessages(library(knitr, warn.conflicts=FALSE))
suppressMessages(library(lmtest, warn.conflicts=FALSE))
suppressMessages(library(MTS, warn.conflicts=FALSE))
suppressMessages(library(quantmod, warn.conflicts=FALSE))
suppressMessages(library(readxl, warn.conflicts=FALSE))
suppressMessages(library(sandwich, warn.conflicts=FALSE))
suppressMessages(library(smooth, warn.conflicts=FALSE))
suppressMessages(library(strucchange, warn.conflicts=FALSE))
suppressMessages(library(timeDate, warn.conflicts=FALSE))
suppressMessages(library(timeSeries, warn.conflicts=FALSE))
suppressMessages(library(trend, warn.conflicts=FALSE))
suppressMessages(library(TSA, warn.conflicts=FALSE))
suppressMessages(library(tscount, warn.conflicts=FALSE))
suppressMessages(library(urca, warn.conflicts=FALSE))
suppressMessages(library(vars, warn.conflicts=FALSE))
suppressMessages(library(zoo, warn.conflicts=FALSE))# Abra o Terminal (Prompt de Comando (CMD) do Windows)
# python.exe -m pip install --upgrade pip
# Digite: pip install pacote
import numpy
import matplotlib
import scipy
import statistics
import pandas
# import math
# import fractions
# import sympy
# import cmathAdicionar
Time Series for Data Science - Woodward et al - 2022
Spatial Predictive Modelling with R - Li - 2022
Wang, T., (2020) “Detection of Learning in Longitudinal Assessment via Score-based Tests”, Practical Assessment, Research, and Evaluation 25(1): 1. doi: https://doi.org/10.7275/gf0x-v588
HH::seqplot: Time series plot,HH::seqplotForecastHH::seqplot(co2)
sunspot.monthwavemulcor
Material
- HTML de R Markdown em
RPubs
Repositório de dados
Geral
Time series - UCI Machine Learning Repository: https://archive.ics.uci.edu/ml/datasets.php?format=&task=&att=&area=&numAtt=&numIns=&type=ts&sort=nameUp&view=table
Saúde
DATASUS:
https://datasus.saude.gov.br/transferencia-de-arquivos/
https://opendatasus.saude.gov.br/dataset/
IBGE - PNS - Pesquisa Nacional de Saúde: https://www.ibge.gov.br/estatisticas/sociais/saude/29540-2013-pesquisa-nacional-de-saude.html?edicao=9163&t=resultados
REDCap: https://redcapbrasil.com.br/
National Center for Health Statistics: https://www.cdc.gov/nchs/data_access/ftp_data.htm
HealthData.gov: https://beta.healthdata.gov/browse
Health Statistics & Data: https://guides.lib.berkeley.edu/publichealth/healthstatistics/usingstats
World Health Organization data collections: https://www.who.int/data/collections
National Cancer Institute: The Surveillance, Epidemiology, and End Results (SEER) Program https://seer.cancer.gov/ https://nccrexplorer.ccdi.cancer.gov/
Sociodemográfico
OECD.Stat: http://stats.oecd.org/
IPEA data: http://www.ipeadata.gov.br
General Social Survey: https://gss.norc.org/
Econômico-financeiro
Finding and using financial data at Princeton University: http://libguides.princeton.edu/c.php?g=84125&p=541948
GoogleFinanças: https://www.google.com/finance
Yahoo!Finanças: http://br.finance.yahoo.com
CBOE DataShop: http://www.marketdataexpress.com/default.aspx
IVolatility DataTools|Data Download Services: http://www.ivolatility.com/data_download.j
Banco Central do Brasil: http://www.bcb.gov.br/htms/selic/selicdiarios.asp
Conteúdo
- Modelo SARIMA
- Função de transferência: SARIMAX
- Análise de intervenção
O que é série temporal?
Série temporal ou sucessão cronológica é uma sequência numérica indexada equispaçadamente no tempo com pelo menos 50 observações.
Cada série é uma variável intervalar com medidas repetidas: delineamento intraparticipantes de sujeito-único.
Objetivo: prever fora do domínio de observação (extrapolação) usando uma combinação linear de valores anteriores.
Análise de série temporal: Saúde e R
Questões
Como simplificar e tornar mais fácil a compreensão de Análise de Série Temporal/ Sucessão Cronológica para a área da Saúde?
Qual a relevância da Análise de Série Temporal para a área da saúde?
Respostas
Para simplificar a análise de dados da área da saúde, decidiu-se pelo uso do R, pois:
- Software livre e gratuito;
- Disponível para LINUX, Windows e MacOSX;
- Muitos pacotes desenvolvidos em R de boa qualidade para analisar dados da área da saúde;
- Expectativa dessa situação continuar por um longo período;
- Scripts R desenvolvidos pelo autor para realizar procedimentos específicos.
Para simplificar a exposição da Análise de Série Temporal, ela é exposta de maneira concisa e são usados exemplos com dados observados para mostrar sua aplicação na área da saúde.
A importância Análise de Série Temporal aumentou substancialmente nos últimos anos, pois:
- Há mais dados secundários da área da saúde disponíveis pela internet numa frequência cada vez maior;
- Hardware e software de Análise de Série Temporal mais poderosos e amigáveis;
- Internet facilita baixar dados e softwares livres e gratuitos.
Econometria Financeira para dados de baixa frequência (diária ou superior)
- A sequência de valores é série temporal, pois os valores são equiespaçados
- Série temporal do retorno segue processo NIID/ GWN
- Analisar principalmente as relações “causais” entre séries temporais (precedência temporal por meio do método de causalidade de Granger) com o uso de cointegração, VAR e VEC
- As séries macroeconômicas raramente passam das poucas centenas de observações (quando, na melhor das hipóteses, se têm observações mensais) e, em geral, a quantidade de dados é insuficiente para os testes econométricos poderosos
- Dados de baixa qualidade
- Séries com periodicidades distintas, em geral
Análise de sujeito único
Morley, 2018, Chapter 5
Morley, S (2018) Single-case methods in clinical psychology: A practical guide. NY: Routledge. Chapter 5: Visual analysis if single-case data; Chapter 6: Statistical analysis of single-case data.
Spata, AV (2005) Métodos de pesquisa: ciências do comportamento e diversidade humana. RJ: LTC. Capítulo 13: Delineamentos de sujeito único.
Velicer, WF & Fava, JL (2003) Time series analysis. In Schinka, J & Velicer, WF (Eds.) Research Methods in Psychology (581-606). Volume 2, Handbook of Psychology (Weiner, IB Editor-inChief). NY: Wiley.
Hamaker, EL (2013) Why researchers should think “within-person”: A paradigmatic rationale. In Mehl MR &, Conner TS (Eds.) (2013) Handbook of research methods for studying daily life. NY: Guilford. p. 43-61.
Delineamento de série temporal
- Unidade observacional: individuada e agregada
- Participação do pesquisador: não intervencional ou intervencional
- Temporalidade: longitudinal com medidas equiespaçadas (periódico)
“A essência do delineamento de série temporal é a presença de um processo periódico de medida em algum grupo ou indivíduo e a introdução de uma mudança experimental nessa série temporal de medidas, cujos resultados são indicador por uma descontinuidade nas medidas registradas na série temporal.”
Campbell & Stanley, 1969, p. 67
Os delineamentos de séries temporais, em que uma mesma área ou população (unidade observacional agregada), é estudada em momentos distintos e periódicos de tempo (e.g., anualmente), são classificados como um subtipo de estudo ecológico (ecológico temporal), sendo que cada unidade de tempo consiste numa unidade ecológica completa.
O estudo com delineamento ecológico tem a capacidade de gerar hipótese temporal.
“Um estudo ecológico centra-se na comparação de grupos, em vez de indivíduos. As variáveis em uma análise ecológica podem ser medidas agregadas, medidas ambientais ou medidas globais. O propósito de uma análise ecológica pode ser fazer inferências biológicas sobre os efeitos nos riscos individuais ou fazer inferências ecológicas sobre os efeitos nas taxas de grupo. Os delineamentos de estudos ecológicos podem ser classificados em duas dimensões: (a) se o grupo primário é medido (estudo exploratório versus analítico); e (b) se os sujeitos são agrupados por local (estudo de grupo múltiplo), por tempo (estudo de tendência temporal) ou por local e tempo (estudo misto). Apesar de várias vantagens práticas dos estudos ecológicos, existem muitos problemas metodológicos que limitam severamente a inferência causal, incluindo viés ecológico e de nível cruzado, problemas de controle de fatores de confusão, classificação incorreta dentro do grupo, falta de dados adequados, ambiguidade temporal, colinearidade e migração entre grupos.”
Morgenstern, 1995, resumo
“Estudo exploratório de série temporal: Um estudo exploratório de tendência temporal ou de série temporal envolve uma comparação das taxas de doença ao longo do tempo em uma população geograficamente definida. Além de fornecer gráficos de tendências temporais, os dados de séries temporais também podem ser usados para prever taxas e tendências futuras. Esta última aplicação, que é mais comum nas ciências sociais do que na epidemiologia, geralmente envolve ajustar os dados de desfecho por modelo de média móvel integrado autorregressivo (ARIMA) (30, 48). O método de modelagem ARIMA também pode ser estendido para avaliar o impacto de uma intervenção populacional (43), para estimar associações entre duas ou mais variáveis de séries temporais (7, 48) e para estimar associações em um delineamento ecológico misto (60; Veja abaixo). Um tipo especial de análise exploratória de tendência temporal frequentemente usada por epidemiologistas é a análise de idade-período-coorte (ou coorte). Por meio de gráficos ou técnicas de modelagem formal, o objetivo desta abordagem é estimar os efeitos separados de três variáveis dependentes do tempo na taxa de doença: idade, período (tempo do calendário) e coorte de nascimento (ano de nascimento) (32, 35). Por causa da dependência linear dessas três variáveis, há uma limitação (problema de identificação) com a interpretação dos resultados da coorte idade-período. O problema é que cada conjunto de dados tem explicações alternativas em relação à combinação de efeitos de idade, período e coorte; não há um conjunto único de parâmetros de efeito quando todas as três variáveis são consideradas simultaneamente. A única maneira de decidir qual interpretação deve ser aceita é considerar os achados à luz do conhecimento prévio e, possivelmente, restringir o modelo ignorando um efeito. Lee et al. (40) conduziram uma análise de idade-período-coorte da mortalidade por melanoma entre homens brancos nos EUA entre 1951 e 1975. Eles concluíram que o aparente aumento na taxa de mortalidade por melanoma era devido principalmente a um efeito de coorte. Ou seja, as pessoas nascidas em anos mais recentes experimentaram ao longo de suas vidas uma taxa mais alta do que as pessoas nascidas mais cedo. Em um artigo subsequente, Lee (39) especulou que esse efeito de coorte pode refletir aumentos na exposição à luz solar ou queimaduras solares durante a juventude.
Estudo analítico de série temporal: Neste tipo de estudo de tendência temporal, avaliamos a associação ecológica entre mudança no nível médio de exposição ou prevalência e mudança na taxa de doença em uma população geograficamente definida. Assim como nos projetos exploratórios, esse tipo de avaliação pode ser feito por simples gráficos ou por modelagem de regressão de séries temporais (por exemplo, 48). Com qualquer uma das abordagens, no entanto, a interpretação dos resultados é muitas vezes complicada por dois problemas. Primeiro, as mudanças na classificação da doença e nos critérios diagnósticos podem produzir resultados muito enganosos. Em segundo lugar, a latência da doença em relação à exposição primária pode ser longa, variável entre os casos ou simplesmente desconhecida. Assim, empregar uma defasagem arbitrária entre as observações - ou uma defasagem empiricamente definida que maximize a associação estimada entre as duas tendências - também pode produzir resultados enganosos (28). Darby & Doll (13) examinaram as associações entre a dose média anual absorvida de precipitação radioativa de testes de armas e a taxa de incidência de leucemia infantil em três países europeus entre 1945 e 1985. Embora a taxa de leucemia tenha variado ao longo do tempo em cada país, eles não encontraram evidências convincentes de que essas mudanças fossem atribuíveis a mudanças na radiação radioativa.”
Morgenstern, 1995, p. 67-8
“Discussão: Quase meio século atrás, Bereiter (1963) escreveu em seu capítulo de abertura do livro Problems in Measuring Change: “Embora seja comum que a pesquisa seja bloqueada por alguma dificuldade na metodologia experimental, realmente não há muitos exemplos nas ciências comportamentais de questões promissoras que não foram pesquisadas devido a deficiências na metodologia estatística. Questões que tratam de mudanças psicológicas podem muito bem constituir as exceções mais importantes”(p. 3). Desde a introdução da modelagem da curva de crescimento latente, esse obstáculo metodológico estatístico foi amplamente superado no que diz respeito aos processos de desenvolvimento. No entanto, só muito recentemente estamos testemunhando o surgimento de soluções para o estudo de processos estacionários. Esta descoberta oportuna e emocionante pode ser parcialmente atribuída a dois desenvolvimentos técnicos recentes. Em primeiro lugar, a introdução de assistentes pessoais digitais, telefones celulares e questionários online facilitou substancialmente a coleta de dados associados à vida diária. Como resultado, mais e mais pesquisadores têm conjuntos de dados que consistem em grandes quantidades de medições repetidas das mesmas pessoas. Em segundo lugar, o maior poder computacional dos computadores permitiu novos e mais complicados modelos de vários níveis pelos quais podemos modelar esses dados longitudinais intensivos. Especificamente, há um novo desenvolvimento que visa a investigação de processos que ocorrem no nível interno, ao mesmo tempo em que toma emprestado a força de outros casos combinando esses modelos internos em novos modelos multiníveis. Com esses desenvolvimentos técnicos em vigor, espera-se que na próxima década veremos um aumento no desenvolvimento e aplicação de novos modelos multiníveis, que são muito mais adequados às hipóteses específicas que os pesquisadores desejam investigar. Agora cabe aos pesquisadores determinar quais questões eles realmente desejam investigar e escolher as ferramentas estatísticas apropriadas para responder a essas questões. Por exemplo, pode-se perguntar não apenas se crianças com mais problemas comportamentais aumentam o nível de rejeição em suas mães, mas também se um aumento nos problemas comportamentais de uma criança leva a um aumento na rejeição materna: Embora a primeira questão seja entre pergunta pessoal, a segunda é uma pergunta interna. Em vez de tentar derivar a resposta para a segunda questão investigando o punho, os pesquisadores agora podem escolher o método apropriado para lidar com sua questão de pesquisa. Os pesquisadores têm feito perguntas pessoais desde o início da psicologia e, por muito tempo, eles não tiveram escolha a não ser fazer a improvável suposição de que os resultados populacionais refletiam fenômenos pessoais. Agora, finalmente, temos acesso a técnicas estatísticas desenvolvidas especificamente para lidar com questões pessoais. Isso implica que a barreira metodológica estatística que Bereiter apontou há quase 50 anos está sendo totalmente superada e estamos no limiar de uma nova e promissora etapa na pesquisa de processos.”
Hamaker, 2013
“A compreensão da dinâmica interacional entre vários processos é um dos desafios mais importantes da psicologia e da medicina psicossomática. Os pesquisadores que exploram o comportamento ou outros fenômenos psicológicos lidam principalmente com dados ordinais ou intervalares. Valores ausentes e medições não equidistantes conseqüentes representam um problema geral de estudos longitudinais deste campo. A maioria das metodologias orientadas a processos foi originalmente projetada para dados equidistantes medidos em escalas de razão. Portanto, o objetivo deste artigo é esclarecer as condições para o desempenho satisfatório dos métodos longitudinais com dados típicos da pesquisa psicológica e psicossomática. Este estudo examina o desempenho do teste de Johansen, um procedimento que incorpora um conjunto de técnicas sofisticadas de séries temporais, em referência à qualidade dos dados utilizando o método de Monte Carlo. Os principais resultados dos estudos de simulação realizados são: (1) As análises de séries temporais requerem amostras de pelo menos 70 observações para uma estimativa e inferência precisas. (2) Dados discretos e falhas na equidistância das medições devido a valores irregulares em falta parecem não ser problemáticos. (3) As características relevantes de processos estacionários podem ser adequadamente capturadas usando escalas ordinais de 5 ou 7 pontos. (4) Para processos de tendências, são necessárias escalas de pelo menos 10 pontos para garantir uma qualidade aceitável de estimativa e inferência.”
Stadnitski, 2020, resumo
“Delineamentos de estudo antes-e-depois são frequentemente usados para quantificar o impacto das intervenções de saúde no nível da população nos processos de cuidado médico e nos desfechos de saúde no nível da população. Eles contam com o “experimento natural” resultante da implementação das intervenções, dividindo o tempo em períodos de “pré-intervenção” e “pós-intervenção”. No entanto, estudos observacionais baseados em um pequeno número de medições pré e pós-intervenção são propensos a viés, pois não levam em conta as tendências subjacentes pré-existentes de curto e longo prazo. Em contraste, a análise de séries temporais interrompidas (ITS) (também chamada de “análise de intervenção”) é mais robusta, pois controla esses problemas rastreando longitudinalmente o resultado antes e depois de uma intervenção. O ITS é considerado um dos melhores desenhos para estabelecer causalidade quando os ensaios clínicos randomizados (ECR) não são viáveis nem éticos. De fato, quando combinados com uma série de controle, os projetos de ITS geralmente geram resultados semelhantes aos ECR. Vários artigos publicados abordaram o tema do uso de abordagens de ITS para avaliar intervenções de saúde. No entanto, eles se concentraram principalmente na regressão segmentada, a forma mais simples de análise de ITS. Modelos de regressão segmentada usam o tempo como preditor variável; um modelo de regressão segmentado simples pode ser expresso como:
\[Y_t = \alpha+\beta_1 \; \text{time} +\beta_2\; \text{intervention} + \beta_3 \;\text{time since intervention} + \epsilon_t\]
Sendo que \(Y_t\) é o resultado em um determinado momento (\(t\)), a variável time representa o tempo desde o início do período de estudo, a variável intervention indica se o tempo \(t\) é antes (0) ou após (1) a implementação da intervenção, e a variável time since intervention representa o tempo decorrido desde a implementação da intervenção, tendo o valor 0 antes da intervenção. Uma suposição chave da regressão linear é que os erros (resíduos) são independentes e não correlacionados. No entanto, esta suposição é frequentemente violada com séries temporais. A abordagem de regressão segmentada é mais apropriada quando uma série temporal tem uma tendência linear ou de outra forma facilmente modelada e resíduos distribuídos independentemente. Na prática, os padrões nos dados podem ser pouco claros ou difíceis de identificar, com variação considerável. Assim, algumas séries temporais podem não ser passíveis de regressão segmentada devido à dificuldade em modelar a estrutura de autocorrelação. Uma alternativa à regressão segmentada são os modelos Autoregressive Integrated Moving Average (ARIMA). Os modelos ARIMA diferem da regressão segmentada, pois o resultado \(Y_t\) é regredido apenas no resultado medido em pontos de tempo anteriores (não no próprio tempo). No entanto, há pouca orientação na literatura sobre como ajustar esses modelos no contexto de análise de ITS. Dada a quantidade e complexidade dos dados de saúde que estão sendo coletados e disponibilizados para pesquisa, o ARIMA tornou-se uma ferramenta cada vez mais útil para pesquisadores interessados em avaliar intervenções em grande escala. Neste artigo, descreveremos a teoria subjacente por trás dos modelos ARIMA e como eles podem ser usados para avaliar intervenções em nível populacional, como a introdução de políticas de saúde, ilustradas usando um exemplo de introdução de uma política de saúde para impedir a prescrição inadequada de quetiapina, um antipsicótico, na Austrália.”
Schaffer et al., 2021, p. 1-2
Livros de análise de série temporal
2015
Livros de referência de análise de série temporal
Livro didático de análise de série temporal
Método de Box-Jenkins
Na prática, raramente encontra-se uma série exatamente conforme um correlograma ACF ou PACF teóricos.
As ambiguidades do método de Box-Jenkins podem levar dois igualmente bons econometristas a estimar e prever uma série usando processos SARIMA bem diferentes.
É prudente evitar o uso de séries já dessazonalizadas.
O gráfico de sequência fornece informações úteis a respeito de outliers, missing values e quebras estruturais.
A ideia fundamental da abordagem de Box-Jenkins é o princípio da parcimônia, pois modelos parcimoniosos preveem melhor que modelos superparametrizados.
Um modelo parcimonioso tem poucos parâmetros e é bem ajustado aos dados.
Quando menor o BIC, mais parcimonioso é o modelo.
Muitos veem a necessidade de julgamento e da experiência do econometrista como uma séria deficiência de um procedimento que é projetado para ser científico.
Conforme Enders (2015), o método de Box_Jenkins é melhor aprendido por meio da experiência, pois análise de série temporal é mais uma arte do que uma ciência.
Airline passengers, 1949 to 1960:
AirPassengers{datasets}
The classic Box, Jenkins & Reinsel (2008) airline data.
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
1952 171 180 193 181 183 218 230 242 209 191 172 194
1953 196 196 236 235 229 243 264 272 237 211 180 201
1954 204 188 235 227 234 264 302 293 259 229 203 229
1955 242 233 267 269 270 315 364 347 312 274 237 278
1956 284 277 317 313 318 374 413 405 355 306 271 306
1957 315 301 356 348 355 422 465 467 404 347 305 336
1958 340 318 362 348 363 435 491 505 404 359 310 337
1959 360 342 406 396 420 472 548 559 463 407 362 405
1960 417 391 419 461 472 535 622 606 508 461 390 432[1] "ts" Time-Series [1:144] from 1949 to 1961: 112 118 132 129 121 135 148 148 136 119 ...plot(1:length(Dados), Dados, type="p", cex=.5, xlab="Mês",
main="Airline passengers, 1949 to 1960",
ylab="Número de passageiros x 10^3", col="black")
lines(1:length(Dados), Dados, type="l", col="gray")
bw <- KernSmooth::dpill(1:length(Dados), Dados, gridsize=length(Dados))
lp <- KernSmooth::locpoly(1:length(Dados), Dados,
bandwidth=bw,
gridsize=length(Dados))
smooth <- lp$y
lines(1:length(Dados), smooth, type="l", lty=3)
Min. 1st Qu. Median Mean 3rd Qu. Max.
104.0 180.0 265.5 280.3 360.5 622.0 $index
[1] 7 19 31 116 128 139 140
$replacements
[1] 185.2686 202.9415 231.6162 452.7610 510.6776 581.2652 549.7057
Created Using changepoint version 2.3
Changepoint type : Change in mean
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 14.90944
Minimum Segment Length : 1
Maximum no. of cpts : Inf
Number of changepoints: 123
NULL
Created Using changepoint version 2.3
Changepoint type : Change in variance
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 14.90944
Minimum Segment Length : 2
Maximum no. of cpts : Inf
Changepoint Locations : 41 101
NULL
Optimal 6-segment partition:
Call:
breakpoints.formula(formula = Dados ~ 1)
Breakpoints at observation number:
26 50 77 98 123
Corresponding to breakdates:
1951(2) 1953(2) 1955(5) 1957(2) 1959(3) 
plot(1:length(Dados), Dados, type="p", cex=.5,
ylim=c(min(Dados)*0.95,max(Dados)*1.05),
main="Airline passengers, 1949 to 1960",
xlab="Mês",
ylab="Número de passageiros x 10^3", col="black")
lines(1:length(Dados), Dados, type="l", col="gray")
for(i in 1:length(sc$breakpoints)) abline(v=sc$breakpoints[i], lty=2)
segments(0, cpm@param.est$mean[1],
cpm@cpts[1], cpm@param.est$mean[1],
col="black")
segments(cpm@cpts[1], cpm@param.est$mean[2],
cpm@cpts[2], cpm@param.est$mean[2],
col="black")
print(cpm@param.est$mean) [1] 112.0000 118.0000 130.5000 121.0000 135.0000 148.0000 136.0000 119.0000
[9] 104.0000 116.5000 126.0000 141.0000 135.0000 125.0000 149.0000 170.0000
[17] 158.0000 133.0000 114.0000 140.0000 147.5000 178.0000 163.0000 172.0000
[25] 178.0000 199.0000 184.0000 162.0000 146.0000 168.5000 180.0000 193.0000
[33] 182.0000 218.0000 230.0000 242.0000 209.0000 191.0000 172.0000 195.3333
[41] 235.5000 229.0000 243.0000 264.0000 272.0000 237.0000 211.0000 180.0000
[49] 202.5000 188.0000 235.0000 227.0000 234.0000 264.0000 302.0000 293.0000
[57] 259.0000 229.0000 203.0000 229.0000 242.0000 233.0000 268.6667 315.0000
[65] 364.0000 347.0000 312.0000 274.0000 237.0000 279.6667 316.0000 374.0000
[73] 413.0000 405.0000 355.0000 306.0000 271.0000 306.0000 315.0000 301.0000
[81] 356.0000 348.0000 355.0000 422.0000 466.0000 404.0000 347.0000 305.0000
[89] 338.0000 318.0000 362.0000 348.0000 363.0000 435.0000 491.0000 505.0000
[97] 404.0000 359.0000 310.0000 337.0000 360.0000 342.0000 406.0000 396.0000
[105] 420.0000 472.0000 548.0000 559.0000 463.0000 407.0000 362.0000 405.0000
[113] 417.0000 391.0000 419.0000 461.0000 472.0000 535.0000 622.0000 606.0000
[121] 508.0000 461.0000 390.0000 432.0000bw <- KernSmooth::dpill(1:length(Dados), Dados, gridsize=length(Dados))
lp <- KernSmooth::locpoly(1:length(Dados), Dados,
bandwidth=bw,
gridsize=length(Dados))
smooth <- lp$y
lines(1:length(Dados), smooth, type="l", lty=1, lwd=1.5)
# https://www.mssqltips.com/sqlservertip/6778/time-series-forecasting-methods-with-r-examples/
fts <- forecast::snaive(Dados, h=60)
summary(fts)
Forecast method: Seasonal naive method
Model Information:
Call: forecast::snaive(y = Dados, h = 60)
Residual sd: 36.3157
Error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 31.77273 36.31574 32.0303 11.12439 11.24871 1 0.7464603
Forecasts:
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
Jan 1961 417 370.4595 463.5405 345.8224 488.1776
Feb 1961 391 344.4595 437.5405 319.8224 462.1776
Mar 1961 419 372.4595 465.5405 347.8224 490.1776
Apr 1961 461 414.4595 507.5405 389.8224 532.1776
May 1961 472 425.4595 518.5405 400.8224 543.1776
Jun 1961 535 488.4595 581.5405 463.8224 606.1776
Jul 1961 622 575.4595 668.5405 550.8224 693.1776
Aug 1961 606 559.4595 652.5405 534.8224 677.1776
Sep 1961 508 461.4595 554.5405 436.8224 579.1776
Oct 1961 461 414.4595 507.5405 389.8224 532.1776
Nov 1961 390 343.4595 436.5405 318.8224 461.1776
Dec 1961 432 385.4595 478.5405 360.8224 503.1776
Jan 1962 417 351.1818 482.8182 316.3397 517.6603
Feb 1962 391 325.1818 456.8182 290.3397 491.6603
Mar 1962 419 353.1818 484.8182 318.3397 519.6603
Apr 1962 461 395.1818 526.8182 360.3397 561.6603
May 1962 472 406.1818 537.8182 371.3397 572.6603
Jun 1962 535 469.1818 600.8182 434.3397 635.6603
Jul 1962 622 556.1818 687.8182 521.3397 722.6603
Aug 1962 606 540.1818 671.8182 505.3397 706.6603
Sep 1962 508 442.1818 573.8182 407.3397 608.6603
Oct 1962 461 395.1818 526.8182 360.3397 561.6603
Nov 1962 390 324.1818 455.8182 289.3397 490.6603
Dec 1962 432 366.1818 497.8182 331.3397 532.6603
Jan 1963 417 336.3895 497.6105 293.7169 540.2831
Feb 1963 391 310.3895 471.6105 267.7169 514.2831
Mar 1963 419 338.3895 499.6105 295.7169 542.2831
Apr 1963 461 380.3895 541.6105 337.7169 584.2831
May 1963 472 391.3895 552.6105 348.7169 595.2831
Jun 1963 535 454.3895 615.6105 411.7169 658.2831
Jul 1963 622 541.3895 702.6105 498.7169 745.2831
Aug 1963 606 525.3895 686.6105 482.7169 729.2831
Sep 1963 508 427.3895 588.6105 384.7169 631.2831
Oct 1963 461 380.3895 541.6105 337.7169 584.2831
Nov 1963 390 309.3895 470.6105 266.7169 513.2831
Dec 1963 432 351.3895 512.6105 308.7169 555.2831
Jan 1964 417 323.9190 510.0810 274.6449 559.3551
Feb 1964 391 297.9190 484.0810 248.6449 533.3551
Mar 1964 419 325.9190 512.0810 276.6449 561.3551
Apr 1964 461 367.9190 554.0810 318.6449 603.3551
May 1964 472 378.9190 565.0810 329.6449 614.3551
Jun 1964 535 441.9190 628.0810 392.6449 677.3551
Jul 1964 622 528.9190 715.0810 479.6449 764.3551
Aug 1964 606 512.9190 699.0810 463.6449 748.3551
Sep 1964 508 414.9190 601.0810 365.6449 650.3551
Oct 1964 461 367.9190 554.0810 318.6449 603.3551
Nov 1964 390 296.9190 483.0810 247.6449 532.3551
Dec 1964 432 338.9190 525.0810 289.6449 574.3551
Jan 1965 417 312.9323 521.0677 257.8422 576.1578
Feb 1965 391 286.9323 495.0677 231.8422 550.1578
Mar 1965 419 314.9323 523.0677 259.8422 578.1578
Apr 1965 461 356.9323 565.0677 301.8422 620.1578
May 1965 472 367.9323 576.0677 312.8422 631.1578
Jun 1965 535 430.9323 639.0677 375.8422 694.1578
Jul 1965 622 517.9323 726.0677 462.8422 781.1578
Aug 1965 606 501.9323 710.0677 446.8422 765.1578
Sep 1965 508 403.9323 612.0677 348.8422 667.1578
Oct 1965 461 356.9323 565.0677 301.8422 620.1578
Nov 1965 390 285.9323 494.0677 230.8422 549.1578
Dec 1965 432 327.9323 536.0677 272.8422 591.1578
Time elapsed: 0.02 seconds
Model estimated using sma() function: SMA(3)
With backcasting initialisation
Distribution assumed in the model: Normal
Loss function type: MSE; Loss function value: 2292.684
ARMA parameters of the model:
Lag 1
AR(1) 0.3333
AR(2) 0.3333
AR(3) 0.3333
Sample size: 144
Number of estimated parameters: 1
Number of degrees of freedom: 143
Information criteria:
AIC AICc BIC BICc
1524.851 1524.879 1527.821 1527.891 Jan Feb Mar Apr May Jun Jul Aug
1961 427.6667 416.5556 425.4074 423.2099 421.7243 423.4472 422.7938 422.6551
Sep Oct
1961 422.9654 422.8047
Box-Pierce test
data: Dados
X-squared = 129.43, df = 1, p-value < 2.2e-16# https://stats.stackexchange.com/questions/104215/difference-between-series-with-drift-and-series-with-trend
print(aTSA::trend.test(Dados)) # H0: ausência de tendência
Approximate Cox-Stuart trend test
data: Dados
D- = 72, p-value < 2.2e-16
alternative hypothesis: data have a decreasing trend
Cox and Stuart Trend test
data: Dados
z = 6.9282, n = 144, p-value = 4.262e-12
alternative hypothesis: monotonic trend
Mann-Kendall trend test
data: Dados
z = 14.382, n = 144, p-value < 2.2e-16
alternative hypothesis: true S is not equal to 0
sample estimates:
S varS tau
8.327000e+03 3.351643e+05 8.098232e-01
Sen's slope
data: Dados
z = 14.382, n = 144, p-value < 2.2e-16
alternative hypothesis: true z is not equal to 0
95 percent confidence interval:
2.294118 2.621622
sample estimates:
Sen's slope
2.451216
Augmented Dickey-Fuller Test
data: Dados
Dickey-Fuller = -7.3186, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary
Approximate Cox-Stuart trend test
data: diff(log(Dados))
D+ = 37, p-value = 0.3609
alternative hypothesis: data have a increasing trend
Cox and Stuart Trend test
data: diff(log(Dados))
z = 0.91733, n = 143, p-value = 0.359
alternative hypothesis: monotonic trend
Mann-Kendall trend test
data: diff(log(Dados))
z = -0.30718, n = 143, p-value = 0.7587
alternative hypothesis: true S is not equal to 0
sample estimates:
S varS tau
-1.770000e+02 3.282717e+05 -1.743842e-02
Sen's slope
data: diff(log(Dados))
z = -0.30718, n = 143, p-value = 0.7587
alternative hypothesis: true z is not equal to 0
95 percent confidence interval:
-0.0005417372 0.0003785882
sample estimates:
Sen's slope
-7.02176e-05
Box-Pierce test
data: diff(log(Dados))
X-squared = 5.7058, df = 1, p-value = 0.01691
Augmented Dickey-Fuller Test
data: diff(log(Dados))
Dickey-Fuller = -6.4313, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary
plot(Dados, col="gray", main="Dados brutos e dessazonalizados")
ds <- forecast::seasadj(decompose(Dados,"multiplicative"))
lines(ds)
Approximate Cox-Stuart trend test
data: diff(log(ds))
D+ = 37, p-value = 0.3609
alternative hypothesis: data have a increasing trend
Cox and Stuart Trend test
data: diff(log(ds))
z = 1.4967, n = 143, p-value = 0.1345
alternative hypothesis: monotonic trend
Mann-Kendall trend test
data: diff(log(ds))
z = -0.66847, n = 143, p-value = 0.5038
alternative hypothesis: true S is not equal to 0
sample estimates:
S varS tau
-3.840000e+02 3.282767e+05 -3.782692e-02
Sen's slope
data: diff(log(ds))
z = -0.66847, n = 143, p-value = 0.5038
alternative hypothesis: true z is not equal to 0
95 percent confidence interval:
-2.012801e-04 9.798408e-05
sample estimates:
Sen's slope
-4.915741e-05
Box-Pierce test
data: diff(log(ds))
X-squared = 7.4856, df = 1, p-value = 0.006219
Augmented Dickey-Fuller Test
data: diff(log(ds))
Dickey-Fuller = -5.9651, Lag order = 5, p-value = 0.01
alternative hypothesis: stationaryBronquiolite
Cunha, GZ (2016) Modelos de séries temporais para dados de contagem. Dissertação de mestrado do PBE-UEM.
COVID-19
Modelagem Matemática do COVID-19: Atualização de 19.06.2020: https://wp.ufpel.edu.br/fentransporte/2020/06/19/modelagem-matematica-do-covid-19-atualizacao-de-19-06-2020/
Tabagismo: Padrão
Velicer & Fava (2003)
# A tibble: 20 × 2
Dia ncigROD
<dbl> <dbl>
1 1 6
2 2 10
3 3 4
4 4 13
5 5 4
6 6 11
7 7 4
8 8 6
9 9 4
10 10 15
11 11 5
12 12 14
13 13 5
14 14 13
15 15 5
16 16 10
17 17 3
18 18 14
19 19 3
20 20 16[1] "tbl_df" "tbl" "data.frame"tibble [20 × 2] (S3: tbl_df/tbl/data.frame)
$ Dia : num [1:20] 1 2 3 4 5 6 7 8 9 10 ...
$ ncigROD: num [1:20] 6 10 4 13 4 11 4 6 4 15 ...
Time Series:
Start = 1
End = 20
Frequency = 1
ncigROD
[1,] 6
[2,] 10
[3,] 4
[4,] 13
[5,] 4
[6,] 11
[7,] 4
[8,] 6
[9,] 4
[10,] 15
[11,] 5
[12,] 14
[13,] 5
[14,] 13
[15,] 5
[16,] 10
[17,] 3
[18,] 14
[19,] 3
[20,] 16[1] "ts" Time-Series [1:20, 1] from 1 to 20: 6 10 4 13 4 11 4 6 4 15 ...
- attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr "ncigROD"plot(1:length(Dados), Dados, type="p", cex=1, xlab="Dia",
ylab="Numero de cigarros", col="black")
lines(1:length(Dados), Dados, type="l", col="gray")
bw <- KernSmooth::dpill(1:length(Dados), Dados, gridsize=length(Dados))
lp <- KernSmooth::locpoly(1:length(Dados), Dados,
bandwidth=bw,
gridsize=length(Dados))
smooth <- lp$y
lines(1:length(Dados), smooth, type="l", lty=3)
# lines(ksmooth(1:length(Dados), Dados,
# kernel="normal",
# bandwidth=bw),
# type='l', col='black', lty=2)
forecast::tsoutliers(Dados)$index
integer(0)
$replacements
numeric(0)Tabagismo: Correlograma
# A tibble: 20 × 2
Dia ncigROD
<dbl> <dbl>
1 1 6
2 2 10
3 3 4
4 4 13
5 5 4
6 6 11
7 7 4
8 8 6
9 9 4
10 10 15
11 11 5
12 12 14
13 13 5
14 14 13
15 15 5
16 16 10
17 17 3
18 18 14
19 19 3
20 20 16plot(Dados$Dia, Dados$ncigROD, ylim=c(0,17),
xlab="Dia", ylab="Numero de cigarros", main="ROD")
lines(Dados$Dia, Dados$ncigROD, type="l", col="gray")
Dados <- cbind(Dados,
dplyr::lag(Dados[,2], n=1),
dplyr::lag(Dados[,2], n=2),
dplyr::lag(Dados[,2], n=3),
dplyr::lag(Dados[,2], n=4),
dplyr::lag(Dados[,2], n=5))
print(Dados) Dia ncigROD ncigROD ncigROD ncigROD ncigROD ncigROD
1 1 6 NA NA NA NA NA
2 2 10 6 NA NA NA NA
3 3 4 10 6 NA NA NA
4 4 13 4 10 6 NA NA
5 5 4 13 4 10 6 NA
6 6 11 4 13 4 10 6
7 7 4 11 4 13 4 10
8 8 6 4 11 4 13 4
9 9 4 6 4 11 4 13
10 10 15 4 6 4 11 4
11 11 5 15 4 6 4 11
12 12 14 5 15 4 6 4
13 13 5 14 5 15 4 6
14 14 13 5 14 5 15 4
15 15 5 13 5 14 5 15
16 16 10 5 13 5 14 5
17 17 3 10 5 13 5 14
18 18 14 3 10 5 13 5
19 19 3 14 3 10 5 13
20 20 16 3 14 3 10 5Dados <- na.omit(Dados)
n <- length(Dados)
# https://www.datacamp.com/tutorial/autocorrelation-r
round(cor(Dados[,2:7])[1,2:6]*(n-1)/n,2)ncigROD.1 ncigROD.2 ncigROD.3 ncigROD.4 ncigROD.5
-0.66 0.66 -0.67 0.57 -0.67 
Velicer & Fava (2003): Figure 23.4
Dados <- readxl::read_excel("Cigarros.xlsx", sheet="ncigarros")
Dados <- subset(Dados, select=c(ncigROD))
Dados <- ts(data=Dados)
acf.est <- acf(Dados, plot=FALSE)
print(acf.est)
Autocorrelations of series 'Dados', by lag
1 2 3 4 5 6 7 8 9 10 11
-0.727 0.681 -0.633 0.502 -0.565 0.546 -0.525 0.555 -0.427 0.391 -0.304
12 13
0.202 -0.302
Partial autocorrelations of series 'Dados', by lag
1 2 3 4 5 6 7 8 9 10 11
-0.727 0.323 -0.149 -0.140 -0.270 0.125 -0.045 0.066 0.193 -0.058 0.159
12 13
-0.136 -0.237 plot(1:(length(acf.est$acf)), acf.est$acf[1:length(acf.est$acf)], type="l", col="gray",
xlab="lag", ylab="Correlacao", main="Correlograma", ylim=c(-1,1))
lines(1:length(pacf.est$acf), pacf.est$acf, type="l", "black")
abline(h=0, lty=2)
abline(v=1, lty=2)
legend("topright", legend=c("ACF","PACF"),
lty=c(1,1), col=c("gray", "black"), bty="n")
Tabagismo: Previsão
# https://www.mssqltips.com/sqlservertip/6778/time-series-forecasting-methods-with-r-examples/
fts <- forecast::naive(Dados, h=7)
print(fts) Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
21 16 5.3343524 26.66565 -0.3117004 32.31170
22 16 0.9164965 31.08350 -7.0682279 39.06823
23 16 -2.4734435 34.47344 -12.2526938 44.25269
24 16 -5.3312952 37.33130 -16.6234007 48.62340
25 16 -7.8491130 39.84911 -20.4740708 52.47407
26 16 -10.1253944 42.12539 -23.9553427 55.95534
27 16 -12.2186511 44.21865 -27.1567026 59.15670
Time elapsed: 0 seconds
Model estimated using sma() function: SMA(3)
With backcasting initialisation
Distribution assumed in the model: Normal
Loss function type: MSE; Loss function value: 31.3287
ARMA parameters of the model:
Lag 1
AR(1) 0.3333
AR(2) 0.3333
AR(3) 0.3333
Sample size: 20
Number of estimated parameters: 1
Number of degrees of freedom: 19
Information criteria:
AIC AICc BIC BICc
127.6482 127.8704 128.6440 128.9768 Time Series:
Start = 21
End = 30
Frequency = 1
[1] 11.00000 10.00000 12.33333 11.11111 11.14815 11.53086 11.26337 11.31413
[9] 11.36946 11.31565
Airline passengers: SARIMA
Dados <- datasets::AirPassengers
# Fit model to first few years of AirPassengers data
aa <- forecast::auto.arima(window(AirPassengers, end=1956+11/12))
summary(aa)Series: window(AirPassengers, end = 1956 + 11/12)
ARIMA(1,1,0)(1,1,0)[12]
Coefficients:
ar1 sar1
-0.2250 -0.2274
s.e. 0.1076 0.1081
sigma^2 = 92.5: log likelihood = -304.98
AIC=615.97 AICc=616.27 BIC=623.22
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.4424643 8.834477 6.351387 0.1376786 2.870884 0.2174955
ACF1
Training set 0.003854251
Box-Pierce test
data: aa$residuals
X-squared = 0.0014261, df = 1, p-value = 0.9699
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.22498 0.10760 -2.0910 0.03653 *
sar1 -0.22737 0.10812 -2.1029 0.03548 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1air.model <- forecast::Arima(window(Dados,
end=1956+11/12),
order=c(0,1,1),
seasonal=list(order=c(0,1,1),
period=12),
lambda=0) # Box-Cox transformation
summary(air.model)Series: window(Dados, end = 1956 + 11/12)
ARIMA(0,1,1)(0,1,1)[12]
Box Cox transformation: lambda= 0
Coefficients:
ma1 sma1
-0.3941 -0.6129
s.e. 0.1173 0.1076
sigma^2 = 0.001556: log likelihood = 148.76
AIC=-291.53 AICc=-291.22 BIC=-284.27
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.3576253 7.89734 5.788344 0.1458472 2.670181 0.1982148 0.05807465
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ma1 -0.39415 0.11726 -3.3613 0.0007759 ***
sma1 -0.61293 0.10758 -5.6973 1.217e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: air.model$residuals
X-squared = 0.041008, df = 1, p-value = 0.8395Modelo SARIMA
Seasonal Autoregressive Integrated Moving Average
Modelo SARIMA da quantidade de passageiros da companhia aérea
Modelo Airline de Newbold & Bos (1990)
SARIMA(p,d,q)(P,D,Q)s
- Autorregressividade de ordem p: p valores anteriores influenciam na predição do valor corrente da série
- Média móvel de ordem q: q desvios da média (resíduos) da série anteriores influenciam na predição do valor corrente da série
- Diferença finita de ordem d
- Período de sazonalidade s
Airline Passengers Model: SARIMA(0,1,1)(0,1,1)12
Dado bruto: \(X_{t}\)
Dado logaritmizado: \(Y_{t} = lnX_{t}\)
Operador de defasagem: \(B^{k}Y_{t}=Y_{t-k}\)
Termo de erro: \(\{u_{t}\} \sim NIID(0,1)\)
\[(1-B)\left(1-B^{12}\right)Y_{t}=\left(1-\theta_{1}B\right) \left(1-\psi_{1}B^{12}\right)u_{t}\]
\[Y_{t} - Y_{t-1} - \left(Y_{t-12}-Y_{t-13}\right)= u_{t}-\theta_{1}u_{t-1}-\psi_{1}u_{t-2}+\theta_{1}\psi_{1}u_{t-13}\]
\[W_{t}=\epsilon_{t}\]
As autocorrelações são calculadas com \(W_{t}\).
O estimador pontual de \(Y\) com \(h\) passos a diante segue o formato de \(W_{t}\).
Modelo SARIMA
How to Configure SARIMA (https://machinelearningmastery.com/sarima-for-time-series-forecasting-in-python/)
Configuring a SARIMA requires selecting hyperparameters for both the trend and seasonal elements of the series.
Trend Elements
There are three trend elements that require configuration.
They are the same as the ARIMA model; specifically:
p: Trend autoregression order.
d: Trend difference order.
q: Trend moving average order.
Seasonal Elements
There are four seasonal elements that are not part of ARIMA that must be configured; they are:
P: Seasonal autoregressive order.
D: Seasonal difference order.
Q: Seasonal moving average order.
s: The number of time steps for a single seasonal period.
Modelo ARIMA
Velicer & Fava (2003)
Velicer & Fava (2003)
- White noise (ruído branco) = ARIMA(0,0,0)
- Random walk (passeio aleatório) = ARIMA(0,1,0)
- Basic exponential smoothing = ARIMA(0,1,1)
- Brown = ARIMA(0,2,2)
- Holt ou double exponential smoothing = ARIMA(0,2,2)
- Damped Holt = ARIMA(0,1,2)
- Damped-Trend = ARIMA(1,1,2)
- Simple Seasonal = ARIMA(0,1,(1,s,s+1))(0,1,0)
- Winters’ Additive = ARIMA(0,1,s+1)(0,1,0)
- Winters’ Multiplicative = não-ARIMA
Análise de função de transferência
A análise de função de transferência generaliza a análise de intervenção, pois uma extensão natural no modelo de intervenção é permitir que a série independente seja diferente de uma variável dicotômica (dummy) determinística.
Na ausência de retroalimentação (feedback), a análise de função de transferência pode ser útil para previsão e teste de hipótese sobre a relação entre duas séries.
Vestuário masculino por catálogo
IBM SPSS Forecasting, 2014
Uma empresa que opera por meio de catálogo, que está interessada em desenvolver um modelo de previsão, coletou dados sobre vendas mensais entre janeiro de 1989 e dezembro de 1998 de vestuário masculino juntamente com algumas séries que podem ser usadas para explicar a variação nas vendas.
Possíveis preditores são a quantidade de catálogos enviados pelo correio, a quantidade de páginas do catálogo, a quantidade de linhas telefônicas abertas para pedidos, o montante gasto em propaganda impressa e quantidade de atendentes do serviço de venda por catálogo.
Esses preditores são úteis para prever as vendas por catálogo de vestuário masculino?
Observação: antes de converter a variável de data, o arquivo deve ser ordenado cronologicamente.
# Regression with ARIMA errors in R
# https://otexts.com/fpp2/regarima.html
# https://robjhyndman.com/acemsforecasting2018/3-Dynamic-Regression.pdf
Dados <- readxl::read_excel("catalog_seasfac_men.xlsx")
Dados <- timeSeries::timeSeries(data=Dados[,2:5],
charvec=timeDate::timeSequence(from="1989-01-01",
to="1998-12-31",
by="month"),
title="Venda de vestuario masculino por
catalogo: 1989-98")
plot(Dados, main="Sales of Men's Clothing: 1989-98")
Pearson's product-moment correlation
data: Dados$phone and Dados$men
t = 14.478, df = 118, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.7245278 0.8563539
sample estimates:
cor
0.799891 
Pearson's product-moment correlation
data: Dados$mail and Dados$men
t = 15.871, df = 118, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.7581174 0.8750341
sample estimates:
cor
0.8252226 
Pearson's product-moment correlation
data: Dados$page and Dados$men
t = 3.1863, df = 118, p-value = 0.001844
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.1076555 0.4385804
sample estimates:
cor
0.2814656 $index
[1] 36 72 96 108 120
$replacements
[1] 13363.19 20699.00 21184.25 22240.57 24289.32$index
integer(0)
$replacements
numeric(0)$index
[1] 12 24 36 48 60 72 84 96 108 120
$replacements
[1] 8446.0 8906.0 9257.5 9010.5 9918.5 10388.0 10230.5 11515.0 11529.0
[10] 11393.0$index
integer(0)
$replacements
numeric(0)# Lambda <- forecast::BoxCox.lambda(Dados)
# ARIMAX model
xreg <- cbind(Dados$mail, Dados$page, Dados$phone)
colnames(xreg) <- c("mail", "page", "phone")
fitauto <- forecast::auto.arima(Dados$men,
stationary=TRUE,
ic="bic",
lambda=1, # Box-Cox transformation
xreg=xreg)
summary(fitauto)Series: Dados$men
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept mail page phone
0.2402 -20635.193 2.3559 32.4338 299.4961
s.e. 0.1027 2501.035 0.2947 21.5280 49.1399
sigma^2 = 9219273: log likelihood = -1129.96
AIC=2271.91 AICc=2272.66 BIC=2288.64
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -1.628825 2972.396 2211.136 -2.599953 13.76932 0.4876605
ACF1
Training set -0.01804527
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 2.4023e-01 1.0273e-01 2.3386 0.01936 *
intercept -2.0635e+04 2.5010e+03 -8.2507 < 2.2e-16 ***
mail 2.3559e+00 2.9468e-01 7.9948 1.298e-15 ***
page 3.2434e+01 2.1528e+01 1.5066 0.13192
phone 2.9950e+02 4.9140e+01 6.0948 1.096e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.039076, df = 1, p-value = 0.8433xreg <- cbind(Dados$mail, Dados$phone)
colnames(xreg) <- c("mail", "phone")
# IBM SPSS 26: Time Series Modeler: Expert Modeler Criteria
fit <- forecast::Arima(Dados$men,
order=c(0,0,0),
seasonal=list(order=c(0,1,0),period=12),
lambda=1, # Box-Cox transformation
xreg=xreg)
summary(fit)Series: Dados$men
Regression with ARIMA(0,0,0)(0,1,0)[12] errors
Box Cox transformation: lambda= 1
Coefficients:
mail phone
1.6166 355.2929
s.e. 0.4770 43.6056
sigma^2 = 11773404: log likelihood = -1031.43
AIC=2068.86 AICc=2069.09 BIC=2076.9
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -143.9287 3224.878 1801.364 -3.091734 11.37345 0.3972864
ACF1
Training set -0.0201345
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
mail 1.61659 0.47704 3.3888 0.000702 ***
phone 355.29290 43.60559 8.1479 3.704e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fit$residuals
X-squared = 0.048648, df = 1, p-value = 0.8254plot(Dados$men, type="l", col="black", lwd=2,
main="diff(Dados$men,12)~
diff(Dados$mail,12)+diff(Dados$phone,12)
Observed=black Fit=blue")
lines(fit$fitted, type="l", col="blue")
Call:
estimatr::lm_robust(formula = diff(Dados$men, 12) ~ diff(Dados$mail,
12) + diff(Dados$phone, 12))
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
(Intercept) -257.483 370.9263 -0.6942 4.891e-01 -992.961 477.996
diff(Dados$mail, 12) 1.829 0.6121 2.9875 3.502e-03 0.615 3.042
diff(Dados$phone, 12) 361.065 35.2824 10.2336 1.832e-17 291.106 431.023
DF
(Intercept) 105
diff(Dados$mail, 12) 105
diff(Dados$phone, 12) 105
Multiple R-squared: 0.4108 , Adjusted R-squared: 0.3996
F-statistic: 52.37 on 2 and 105 DF, p-value: < 2.2e-16fitauto <- forecast::auto.arima(Dados$men,
stationary=TRUE,
ic="bic",
lambda=1, # Box-Cox transformation
xreg=xreg)
summary(fitauto)Series: Dados$men
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept mail phone
0.288 -19116.166 2.4446 304.9904
s.e. 0.095 2396.267 0.2895 49.4224
sigma^2 = 9318809: log likelihood = -1131.13
AIC=2272.27 AICc=2272.8 BIC=2286.21
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -1.687006 3001.363 2219.18 -2.524513 13.81259 0.4894347
ACF1
Training set -0.02575663
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 2.8798e-01 9.5001e-02 3.0314 0.002434 **
intercept -1.9116e+04 2.3963e+03 -7.9775 1.494e-15 ***
mail 2.4446e+00 2.8950e-01 8.4441 < 2.2e-16 ***
phone 3.0499e+02 4.9422e+01 6.1711 6.782e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.079608, df = 1, p-value = 0.7778xreg <- cbind(Dados$phone)
colnames(xreg) <- c("phone")
fitauto <- forecast::auto.arima(Dados$men,
stationary=TRUE,
ic="bic",
lambda=1, # Box-Cox transformation
xreg=xreg)
summary(fitauto)Series: Dados$men
Regression with ARIMA(0,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
intercept phone
-4395.276 599.055
s.e. 1470.126 41.030
sigma^2 = 14595315: log likelihood = -1159.04
AIC=2324.07 AICc=2324.28 BIC=2332.44
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 4.873442e-12 3788.411 2862.633 -4.062216 17.29913 0.6313467
ACF1
Training set 0.07101713
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
intercept -4395.28 1470.13 -2.9897 0.002792 **
phone 599.05 41.03 14.6004 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.60521, df = 1, p-value = 0.4366xreg <- cbind(Dados$mail)
colnames(xreg) <- c("mail")
fitauto <- forecast::auto.arima(Dados$men,
stationary=TRUE,
ic="bic",
lambda=1, # Box-Cox transformation
xreg=xreg)
summary(fitauto)Series: Dados$men
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept mail
0.2589 -19627.281 3.5345
s.e. 0.0885 2513.584 0.2423
sigma^2 = 12165493: log likelihood = -1147.64
AIC=2303.27 AICc=2303.62 BIC=2314.42
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -6.813439 3444.032 2666.117 -4.371289 17.70102 0.5880055
ACF1
Training set -0.009283068
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 2.5886e-01 8.8537e-02 2.9237 0.003459 **
intercept -1.9627e+04 2.5136e+03 -7.8085 5.788e-15 ***
mail 3.5345e+00 2.4232e-01 14.5861 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.010341, df = 1, p-value = 0.919Melanoma & Mancha Solar
“Houghton, Munster e Viola (1978) mostraram que a taxa de incidência de melanoma maligno ajustada à idade no estado de Connecticut aumentou desde 1935 e que se sobrepõe ao aumento são períodos de 3-5 anos em que o aumento na taxa de incidência é excessivo.
Esses períodos têm um ciclo de 8 a 11 anos e seguem períodos de atividade máxima das manchas solares.
A relação entre os ciclos solares e o melanoma sustenta a hipótese de que o melanoma está relacionado à exposição solar e fornece evidências de que a radiação solar pode desencadear o desenvolvimento de melanoma clinicamente aparente.”
Andrews & Herzberger, 1985, p. 199
Andrews & Herzberger, 1985
Houghton et al., 1978
“Este método revela que existem fortes associações numéricas, com correlações > 0.5, entre essas pequenas, mas distintas tendências de longo prazo no ciclo solar e câncer de pele.”
Valachovic & Zurbenko, 2014
Valachovic & Zurbenko, 2014
Valachovic & Zurbenko, 2014
Fonte dos dados:
- Andrews & Herzberger, 1985, 199-201, Incidence of Malignant Melanoma After Peaks of Sunspot Activity
Dados <- readxl::read_excel("Melanoma.xlsx")
Dados <- timeSeries::timeSeries(data=Dados[,2:4],
charvec=timeDate::timeSequence(from="1936",
to="1972",
by="year"),
title="Age-adjusted incidence of malignant
melanoma in the State of Connecticut
per 100 000 persons")
Dados <- cbind(Dados,diff(log(Dados$Male)))
Dados <- cbind(Dados,log(Dados$Sunspot))
colnames(Dados) <- c("Male", "Total", "Sunspot", "MaleDiffLn", "SunspotLn")
print(Dados)GMT
Male Total Sunspot MaleDiffLn SunspotLn
1936-08-08 1.0 0.9 4.0 -0.22314355 1.3862944
1937-08-08 0.8 0.8 11.5 0.00000000 2.4423470
1938-08-08 0.8 0.8 10.0 0.55961579 2.3025851
1939-08-08 1.4 1.3 8.0 -0.15415068 2.0794415
1940-08-08 1.2 1.4 6.0 -0.18232156 1.7917595
1941-08-08 1.0 1.2 4.0 0.40546511 1.3862944
1942-08-08 1.5 1.7 2.3 0.23638878 0.8329091
1943-08-08 1.9 1.8 1.0 -0.23638878 0.0000000
1944-08-08 1.5 1.6 1.0 0.00000000 0.0000000
1945-08-08 1.5 1.5 2.5 0.00000000 0.9162907
1946-08-08 1.5 1.5 7.5 0.06453852 2.0149030
1947-08-08 1.6 2.0 14.5 0.11778304 2.6741486
1948-08-08 1.8 2.5 13.0 0.44183275 2.5649494
1949-08-08 2.8 2.7 13.0 -0.11332869 2.5649494
1950-08-08 2.5 2.9 8.0 0.00000000 2.0794415
1951-08-08 2.5 2.5 6.5 -0.04082199 1.8718022
1952-08-08 2.4 3.1 2.0 -0.13353139 0.6931472
1953-08-08 2.1 2.4 1.0 -0.10008346 0.0000000
1954-08-08 1.9 2.2 0.5 0.23361485 -0.6931472
1955-08-08 2.4 2.9 1.0 0.00000000 0.0000000
1956-08-08 2.4 2.5 6.0 0.08004271 1.7917595
1957-08-08 2.6 2.6 19.0 0.00000000 2.9444390
1958-08-08 2.6 3.2 18.0 0.52609310 2.8903718
1959-08-08 4.4 3.8 17.5 -0.04652002 2.8622009
1960-08-08 4.2 4.2 12.0 -0.10008346 2.4849066
1961-08-08 3.8 3.9 5.0 -0.11122564 1.6094379
1962-08-08 3.4 3.7 3.5 0.05715841 1.2527630
1963-08-08 3.6 3.3 2.0 0.13005313 0.6931472
1964-08-08 4.1 3.7 1.0 -0.10265415 0.0000000
1965-08-08 3.7 3.9 1.5 0.12675171 0.4054651
1966-08-08 4.2 4.1 3.0 -0.02409755 1.0986123
1967-08-08 4.1 3.8 6.0 0.00000000 1.7917595
1968-08-08 4.1 4.7 10.5 -0.02469261 2.3513753
1969-08-08 4.0 4.4 10.5 0.26236426 2.3513753
1970-08-08 5.2 4.8 10.5 0.01904819 2.3513753
1971-08-08 5.3 4.8 8.0 0.00000000 2.0794415
1972-08-08 5.3 4.8 6.5 -0.22314355 1.8718022
plot(Dados$MaleDiffLn, type="l", ylab="")
lines(Dados$SunspotLn/10, type="l", lty=2, lwd=2)
legend("topright", legend=c("MenDiffLn","SunspotLn"), lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$MaleDiffLn and Dados$SunspotLn
t = 1.0735, df = 35, p-value = 0.2904
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1544063 0.4750800
sample estimates:
cor
0.1785454 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Male,
stationary=TRUE,
ic="bic",
lambda=0, # log: Box-Cox transformation
xreg=Dados$Sunspot)
summary(fitauto)Series: Dados$Male
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 0
Coefficients:
ar1 intercept xreg
0.9535 0.8506 0.0001
s.e. 0.0465 0.4853 0.0088
sigma^2 = 0.04263: log likelihood = 6.23
AIC=-4.47 AICc=-3.22 BIC=1.97
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.1369491 0.4842972 0.3191064 2.005486 13.15987 1.016622
ACF1
Training set -0.09344412
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 9.5354e-01 4.6465e-02 20.5218 < 2e-16 ***
intercept 8.5060e-01 4.8527e-01 1.7529 0.07963 .
xreg 6.2899e-05 8.8479e-03 0.0071 0.99433
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.22369, df = 1, p-value = 0.6362Fonte dos dados:
National Cancer Institute: https://nccrexplorer.ccdi.cancer.gov/
- Melanoma of the Skin
- Recent Trends in SEER Age-Adjusted Incidence Rates, 2004-2019
- By Sex, Observed SEER Incidence Rate, Non-Hispanic White, Ages 40-64, Localized
Solar Influences Data Analysis Center: http://www.sidc.oma.be/
Dados <- readxl::read_excel("Melanoma_Sex.xlsx", sheet="MelanoMale")
Dados <- timeSeries::timeSeries(data=Dados[,2:3],
charvec=timeDate::timeSequence(from="1999",
to="2018",
by="year"),
title="Age-adjusted incidence of malignant
melanoma in the USA: Male")
Dados <- cbind(Dados,log(Dados$Sunspot))
Dados <- cbind(Dados,diff(log(Dados$Melanoma.Male)))
colnames(Dados) <- c("Sunspot", "Melanoma.Male", "SunspotLn", "Melanoma.MaleDiffLn")
print(Dados)GMT
Sunspot Melanoma.Male SunspotLn Melanoma.MaleDiffLn
1999-08-08 136.3 113.8 4.914858 0.104217508
2000-08-08 173.9 126.3 5.158480 0.034244591
2001-08-08 170.4 130.7 5.138149 -0.030293579
2002-08-08 163.6 126.8 5.097424 -0.070232937
2003-08-08 99.3 118.2 4.598146 0.210543351
2004-08-08 65.3 145.9 4.178992 0.028380283
2005-08-08 45.8 150.1 3.824284 -0.117200261
2006-08-08 24.7 133.5 3.206803 0.035323761
2007-08-08 12.6 138.3 2.533697 0.002166848
2008-08-08 4.2 138.6 1.435085 0.008620743
2009-08-08 4.8 139.8 1.568616 0.046129772
2010-08-08 24.9 146.4 3.214868 -0.085522173
2011-08-08 80.8 134.4 4.391977 0.042249547
2012-08-08 84.5 140.2 4.436752 -0.061025915
2013-08-08 94.0 131.9 4.543295 0.020263357
2014-08-08 113.3 134.6 4.730039 0.079928402
2015-08-08 69.8 145.8 4.245634 -0.037740328
2016-08-08 39.8 140.4 3.683867 -0.017241806
2017-08-08 21.7 138.0 3.077312 -0.106165993
2018-08-08 7.0 124.1 1.945910 0.104217508
plot(Dados$Melanoma.MaleDiffLn, type="l", ylab="")
lines(Dados$Sunspot/1000, type="l", lty=2, lwd=2)
legend("topright", legend=c("Melanoma.MaleDiffLn","SunspotLn"),
lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$Melanoma.MaleDiffLn and Dados$SunspotLn
t = 0.12972, df = 18, p-value = 0.8982
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4176071 0.4667696
sample estimates:
cor
0.0305614 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Melanoma.Male,
stationary=TRUE,
ic="bic",
lambda=0, # log: Box-Cox transformation
xreg=Dados$Sunspot)
summary(fitauto)Series: Dados$Melanoma.Male
Regression with ARIMA(0,0,0) errors
Box Cox transformation: lambda= 0
Coefficients:
intercept xreg
4.9485 -6e-04
s.e. 0.0222 3e-04
sigma^2 = 0.004099: log likelihood = 27.64
AIC=-49.29 AICc=-47.79 BIC=-46.3
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.2438441 8.026791 5.912802 -0.1873274 4.470972 0.736677
ACF1
Training set 0.00358662
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
intercept 4.94853643 0.02224107 222.4954 < 2e-16 ***
xreg -0.00064768 0.00033182 -1.9519 0.05095 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.00030616, df = 1, p-value = 0.986Dados <- readxl::read_excel("Melanoma_Sex.xlsx", sheet="MelanoFemale")
Dados <- timeSeries::timeSeries(data=Dados[,2:3],
charvec=timeDate::timeSequence(from="1999",
to="2018",
by="year"),
title="Age-adjusted incidence of malignant
melanoma in the USA: Female")
Dados <- cbind(Dados,log(Dados$Sunspot))
Dados <- cbind(Dados,diff(log(Dados$Melanoma.Female)))
colnames(Dados) <- c("Sunspot", "Melanoma.Female", "SunspotLn", "Melanoma.FemaleDiffLn")
print(Dados)GMT
Sunspot Melanoma.Female SunspotLn Melanoma.FemaleDiffLn
1999-08-08 136.3 153.8 4.914858 0.107763348
2000-08-08 173.9 171.3 5.158480 0.068253154
2001-08-08 170.4 183.4 5.138149 -0.016492952
2002-08-08 163.6 180.4 5.097424 -0.048263695
2003-08-08 99.3 171.9 4.598146 0.084195704
2004-08-08 65.3 187.0 4.178992 0.199989861
2005-08-08 45.8 228.4 3.824284 -0.035200783
2006-08-08 24.7 220.5 3.206803 -0.037897326
2007-08-08 12.6 212.3 2.533697 0.002352389
2008-08-08 4.2 212.8 1.435085 0.038714512
2009-08-08 4.8 221.2 1.568616 -0.034494107
2010-08-08 24.9 213.7 3.214868 -0.049890442
2011-08-08 80.8 203.3 4.391977 0.049422386
2012-08-08 84.5 213.6 4.436752 0.017173780
2013-08-08 94.0 217.3 4.543295 0.024097552
2014-08-08 113.3 222.6 4.730039 0.052931664
2015-08-08 69.8 234.7 4.245634 -0.028523974
2016-08-08 39.8 228.1 3.683867 -0.007922577
2017-08-08 21.7 226.3 3.077312 -0.006205694
2018-08-08 7.0 224.9 1.945910 0.107763348
plot(Dados$Melanoma.FemaleDiffLn, type="l", ylab="")
lines(Dados$Sunspot/1000, type="l", lty=2, lwd=2)
legend("topright", legend=c("Melanoma.FemaleDiffLn","SunspotLn"),
lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$Melanoma.FemaleDiffLn and Dados$SunspotLn
t = 0.60496, df = 18, p-value = 0.5528
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3214369 0.5493664
sample estimates:
cor
0.1411633 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Melanoma.Female,
stationary=TRUE,
ic="bic",
lambda=0, # log: Box-Cox transformation
xreg=Dados$Sunspot)
summary(fitauto)Series: Dados$Melanoma.Female
Regression with ARIMA(0,0,1) errors
Box Cox transformation: lambda= 0
Coefficients:
ma1 intercept xreg
0.9363 5.4118 -0.0013
s.e. 0.2426 0.0387 0.0005
sigma^2 = 0.004114: log likelihood = 27.17
AIC=-46.33 AICc=-43.67 BIC=-42.35
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 1.15578 11.38972 9.205783 0.1797207 4.652125 0.96263 -0.01202982
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ma1 0.93632329 0.24259865 3.8596 0.0001136 ***
intercept 5.41180298 0.03872320 139.7561 < 2.2e-16 ***
xreg -0.00133044 0.00046969 -2.8326 0.0046176 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.074745, df = 1, p-value = 0.7845Dados <- readxl::read_excel("Melanoma_Sex.xlsx", sheet="MelanoMaleSEER")
Dados <- timeSeries::timeSeries(data=Dados[,2:3],
charvec=timeDate::timeSequence(from="2004",
to="2019",
by="year"),
title="Age-adjusted incidence of malignant
melanoma in the USA: SEER Male")
Dados <- cbind(Dados,log(Dados$Sunspot))
Dados <- cbind(Dados,diff(log(Dados$Melanoma.Male)))
colnames(Dados) <- c("Sunspot", "Melanoma.Male", "SunspotLn", "Melanoma.MaleDiffLn")
print(Dados)GMT
Sunspot Melanoma.Male SunspotLn Melanoma.MaleDiffLn
2004-08-08 65.3 31.7 4.178992 0.058209386
2005-08-08 45.8 33.6 3.824284 -0.014992785
2006-08-08 24.7 33.1 3.206803 -0.018293193
2007-08-08 12.6 32.5 2.533697 0.021309787
2008-08-08 4.2 33.2 1.435085 -0.024391453
2009-08-08 4.8 32.4 1.568616 -0.060431739
2010-08-08 24.9 30.5 3.214868 0.013029500
2011-08-08 80.8 30.9 4.391977 0.044311046
2012-08-08 84.5 32.3 4.436752 -0.031449133
2013-08-08 94.0 31.3 4.543295 0.079796917
2014-08-08 113.3 33.9 4.730039 0.000000000
2015-08-08 69.8 33.9 4.245634 0.017544310
2016-08-08 39.8 34.5 3.683867 0.002894358
2017-08-08 21.7 34.6 3.077312 -0.044320400
2018-08-08 7.0 33.1 1.945910 -0.003025721
2019-08-08 3.6 33.0 1.280934 0.058209386
plot(Dados$Melanoma.MaleDiffLn, type="l", ylab="")
lines(Dados$Sunspot/1000, type="l", lty=2, lwd=2)
legend("topright", legend=c("Melanoma.MaleDiffLn","SunspotLn"),
lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$Melanoma.MaleDiffLn and Dados$SunspotLn
t = 1.1545, df = 14, p-value = 0.2676
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2352477 0.6897376
sample estimates:
cor
0.2948399 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Melanoma.Male,
stationary=TRUE,
ic="bic",
lambda=0, # log: Box-Cox transformation
xreg=Dados$Sunspot)
summary(fitauto)Series: Dados$Melanoma.Male
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 0
Coefficients:
ar1 intercept xreg
0.4742 3.4924 -1e-04
s.e. 0.2150 0.0198 4e-04
sigma^2 = 0.001208: log likelihood = 32.58
AIC=-57.16 AICc=-53.53 BIC=-54.07
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set 0.05094932 1.018855 0.8470739 0.05936949 2.593228 0.9011425
ACF1
Training set -0.0004074992
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 0.47420779 0.21504359 2.2052 0.02744 *
intercept 3.49236048 0.01979330 176.4416 < 2e-16 ***
xreg -0.00010457 0.00038832 -0.2693 0.78771
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.00017816, df = 1, p-value = 0.9894Dados <- readxl::read_excel("Melanoma_Sex.xlsx", sheet="MelanoFemaleSEER")
Dados <- timeSeries::timeSeries(data=Dados[,2:3],
charvec=timeDate::timeSequence(from="2004",
to="2019",
by="year"),
title="Age-adjusted incidence of malignant
melanoma in the USA: SEER Female")
Dados <- cbind(Dados,log(Dados$Sunspot))
Dados <- cbind(Dados,diff(log(Dados$Melanoma.Female)))
colnames(Dados) <- c("Sunspot", "Melanoma.Female", "SunspotLn", "Melanoma.FemaleDiffLn")
print(Dados)GMT
Sunspot Melanoma.Female SunspotLn Melanoma.FemaleDiffLn
2004-08-08 65.3 25.2 4.178992 0.057819571
2005-08-08 45.8 26.7 3.824284 0.000000000
2006-08-08 24.7 26.7 3.206803 0.033152207
2007-08-08 12.6 27.6 2.533697 -0.021978907
2008-08-08 4.2 27.0 1.435085 0.036367644
2009-08-08 4.8 28.0 1.568616 -0.028987537
2010-08-08 24.9 27.2 3.214868 -0.007380107
2011-08-08 80.8 27.0 4.391977 0.071458964
2012-08-08 84.5 29.0 4.436752 -0.038669141
2013-08-08 94.0 27.9 4.543295 0.085815920
2014-08-08 113.3 30.4 4.730039 0.060624622
2015-08-08 69.8 32.3 4.245634 0.015361285
2016-08-08 39.8 32.8 3.683867 0.006079046
2017-08-08 21.7 33.0 3.077312 -0.068992871
2018-08-08 7.0 30.8 1.945910 0.084030749
2019-08-08 3.6 33.5 1.280934 0.057819571
plot(Dados$Melanoma.FemaleDiffLn, type="l", ylab="")
lines(Dados$Sunspot/1000, type="l", lty=2, lwd=2)
legend("topright", legend=c("Melanoma.FemaleDiffLn","SunspotLn"),
lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$Melanoma.FemaleDiffLn and Dados$SunspotLn
t = 0.39427, df = 14, p-value = 0.6993
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4123318 0.5708447
sample estimates:
cor
0.1047926 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Melanoma.Female,
stationary=TRUE,
ic="bic",
lambda=0, # log: Box-Cox transformation
xreg=Dados$Sunspot)
summary(fitauto)Series: Dados$Melanoma.Female
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 0
Coefficients:
ar1 intercept xreg
0.8820 3.3865 -5e-04
s.e. 0.1177 0.0772 6e-04
sigma^2 = 0.002814: log likelihood = 25.19
AIC=-42.38 AICc=-38.74 BIC=-39.29
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.3525914 1.41372 1.173911 1.014339 3.982713 0.9728545 -0.3395354
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 0.88197101 0.11773651 7.4911 6.832e-14 ***
intercept 3.38650143 0.07721872 43.8560 < 2.2e-16 ***
xreg -0.00048707 0.00058289 -0.8356 0.4034
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 1.9703, df = 1, p-value = 0.1604Câncer pulmonar & Mancha Solar
“Câncer pulmonar.
Os pacientes com câncer pulmonar registrados investigados totalizaram 1.133 pessoas com 969 homens e 164 mulheres, incluindo câncer pulmonar em diferentes posições como pulmão esquerdo, pulmão direito, folha superior, folha inferior, centro do pulmão, periferia do pulmão e brônquios e câncer pulmonar de diferentes tipos, como adenocarcinoma de pulmão, carcinoma de células alveolares de pulmão e carcinoma de células escamosas.
Comparar a curva de distribuição dos anos de nascimento para os pacientes com câncer de pulmão com a curva de variação periódica da atividade solar, eles têm a mesma semelhança com o câncer gástrico e hepático.
A curva de distribuição do ano de nascimento nos anos de 1926 a 1950 para os pacientes com câncer de pulmão tem uma inter-relação óbvia com a curva de variação periódica da atividade solar (Fig. 4).
Os picos da curva de distribuição do ano de nascimento aparecem nos segmentos de pico e vale da curva de variação periódica da atividade solar, enquanto o número de pacientes com câncer de pulmão, que nasceram nos anos com atividade solar moderada como sunpots entre cerca de 20 a 70, também diminui significativamente. Pelo tratamento estatístico do teste t, a diferença entre eles é significativa com p<0.05.
Somente entre os anos de 1939 e 1942, a curva decrescente do número de pacientes com câncer de pulmão é muito acentuada, e o pico entre 1942 e 1944 é muito baixo.
Possivelmente, isso se deve ao fato de que as pessoas nascidas nestes anos tinham em torno de 50 anos de idade, que se encontra em um estágio de transição de baixa para alta morbidade para o câncer de pulmão, e o número de pacientes é mais afetado pelo fator idade.”
Valachovic & Zurbenko, 2014
Zhang et al., 1993
Fonte dos dados:
National Cancer Institute: https://nccrexplorer.ccdi.cancer.gov/
- Adenocarcinoma of the Lung and Bronchus
- Long-Term Trends in SEER Age-Adjusted Incidence Rates, 1975-2019
- Observed SEER Incidence Rate By Sex, White (includes Hispanic), Ages 40-64
Solar Influences Data Analysis Center: http://www.sidc.oma.be/
Dados <- readxl::read_excel("Lung_Sex.xlsx", sheet="LungMen")
Dados <- timeSeries::timeSeries(data=Dados[,2:3],
charvec=timeDate::timeSequence(from="1975",
to="2019",
by="year"),
title="Age-adjusted incidence of Lung cancer
in USA")
Dados <- cbind(Dados,diff(log(Dados$Lung)))
Dados <- cbind(Dados,diff(log(Dados$Sunspot)))
colnames(Dados) <- c("Sunspot", "Lung", "LungDiffLn", "SunspotDiffLn")
print(Dados)GMT
Sunspot Lung LungDiffLn SunspotDiffLn
1975-08-08 22.5 22.0 0.170892861 -0.201164645
1976-08-08 18.4 26.1 0.080926490 0.758873854
1977-08-08 39.3 28.3 0.110295316 1.203972804
1978-08-08 131.0 31.6 -0.068766857 0.518884665
1979-08-08 220.1 29.5 0.016807118 -0.005466984
1980-08-08 218.9 30.0 0.055119299 -0.095812819
1981-08-08 198.9 31.7 -0.106412594 -0.202739758
1982-08-08 162.4 28.5 -0.024868067 -0.579202921
1983-08-08 91.0 27.8 0.066111025 -0.408216141
1984-08-08 60.5 29.7 0.080819407 -1.077352289
1985-08-08 20.6 32.2 -0.060818740 -0.330663895
1986-08-08 14.8 30.3 0.063919518 0.828787834
1987-08-08 33.9 32.3 0.006172859 1.288769341
1988-08-08 123.0 32.5 -0.063513406 0.540147599
1989-08-08 211.1 30.5 0.022691411 -0.095878792
1990-08-08 191.8 31.2 -0.025975486 0.058229558
1991-08-08 203.3 30.4 -0.023295563 -0.424333592
1992-08-08 133.0 29.7 -0.099020759 -0.558300863
1993-08-08 76.1 26.9 0.000000000 -0.527610470
1994-08-08 44.9 26.9 -0.061321891 -0.581569949
1995-08-08 25.1 25.3 -0.028057953 -0.771862748
1996-08-08 11.6 24.6 -0.008163311 0.912836497
1997-08-08 28.9 24.4 0.036221263 1.116898512
1998-08-08 88.3 25.3 -0.086652117 0.434118231
1999-08-08 136.3 23.2 -0.076099344 0.243622083
2000-08-08 173.9 21.5 0.000000000 -0.020331807
2001-08-08 170.4 21.5 -0.009345862 -0.040724190
2002-08-08 163.6 21.3 -0.130183549 -0.499278853
2003-08-08 99.3 18.7 -0.119120828 -0.419153535
2004-08-08 65.3 16.6 0.000000000 -0.354707945
2005-08-08 45.8 16.6 0.011976191 -0.617480847
2006-08-08 24.7 16.8 -0.018018506 -0.673106430
2007-08-08 12.6 16.5 0.018018506 -1.098612289
2008-08-08 4.2 16.8 0.040821995 0.133531393
2009-08-08 4.8 17.5 -0.040821995 1.646251886
2010-08-08 24.9 16.8 0.133531393 1.177109162
2011-08-08 80.8 19.2 -0.064538521 0.044774569
2012-08-08 84.5 18.0 0.011049836 0.106543248
2013-08-08 94.0 18.2 -0.068208250 0.186744386
2014-08-08 113.3 17.0 0.079137321 -0.484405158
2015-08-08 69.8 18.4 -0.115069330 -0.561767097
2016-08-08 39.8 16.4 -0.024692613 -0.606554652
2017-08-08 21.7 16.0 -0.098440073 -1.131402111
2018-08-08 7.0 14.5 -0.164549387 -0.664976304
2019-08-08 3.6 12.3 0.170892861 -0.201164645
plot(Dados$LungDiffLn, type="l", ylab="")
lines(Dados$SunspotDiffLn/10, type="l", lty=2, lwd=2)
legend("top", legend=c("LungDiffLn","SunspotDiffLn"), lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$LungDiffLn and Dados$SunspotDiffLn
t = 1.6752, df = 43, p-value = 0.1012
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.04962391 0.50440237
sample estimates:
cor
0.2475154 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Lung,
stationary=TRUE,
ic="bic",
lambda=0, # Box-Cox transformation
xreg=Dados$Sunspot)
summary(fitauto)Series: Dados$Lung
Regression with ARIMA(0,0,3) errors
Box Cox transformation: lambda= 0
Coefficients:
ma1 ma2 ma3 intercept xreg
0.9883 0.9923 0.8214 3.0085 0.0011
s.e. 0.1274 0.1923 0.2296 0.0708 0.0006
sigma^2 = 0.01104: log likelihood = 38.08
AIC=-64.15 AICc=-61.94 BIC=-53.31
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.3418257 2.324902 1.877017 -0.3196853 8.136877 1.374189 0.2409149
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ma1 0.98833063 0.12739306 7.7581 8.620e-15 ***
ma2 0.99233236 0.19232575 5.1596 2.474e-07 ***
ma3 0.82140480 0.22955843 3.5782 0.000346 ***
intercept 3.00848500 0.07076109 42.5161 < 2.2e-16 ***
xreg 0.00108440 0.00055485 1.9544 0.050655 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 2.7469, df = 1, p-value = 0.09744Câncer pulmonar & Tabagismo
“Ao discutir suas descobertas, Doll e Hill raciocinaram cuidadosamente ao pesar as evidências a favor e contra o tabagismo como “um fator importante na causa do carcinoma do pulmão”. ”
Schlesselman, 2006, p. 254
“Calculamos os riscos de morte por tabagismo a longo prazo para indivíduos de várias idades e status de tabagismo em termos do excesso de mortalidade contribuído pelo tabagismo, além da mortalidade basal pelas mesmas doenças causadas por outros fatores que não o tabagismo, usando procedimentos de tabela de vida padrão.
Como os dados de mortalidade para categorias específicas de fumantes estavam disponíveis apenas em estudos prospectivos no final da década de 1950, nós os escalamos para os níveis de mortalidade de 1982.
Assumimos, para câncer de pulmão, que as taxas de mortalidade para não fumantes não mudaram e, para outras doenças relacionadas ao tabagismo, que os riscos de morte para fumantes em relação aos não fumantes não mudaram desde a década de 1950.
Probabilidades que resultam de hipóteses alternativas também foram investigadas e são apresentadas.
Cerca de um terço dos fumantes pesados de 35 anos morrerão antes dos 85 anos de doenças causadas pelo tabagismo.
As probabilidades de morte por tabagismo quando comparadas com outras causas podem ser persuasivas como ferramentas de educação pública.
Seu uso efetivo para esse fim é afetado não apenas pelas deficiências no conhecimento factual do público sobre a magnitude dos riscos do tabagismo, mas também por vários aparentes equívocos relacionados à interpretação das informações de risco.
Os dados de risco devem ser apresentados ao público de uma maneira que esclareça esses equívocos e facilite sua compreensão do risco esmagador imposto pelo tabagismo.”
Mattson et al., 1987, resumo
Hecht, 1999
Schlesselman, 2006
Mattson et al., 1987
Mattson et al., 1987
Fonte dos dados:
National Cancer Institute: https://nccrexplorer.ccdi.cancer.gov/
- Adenocarcinoma of the Lung and Bronchus
- Long-Term Trends in SEER Age-Adjusted Incidence Rates, 1975-2019
- Observed SEER Incidence Rate By Sex, White (includes Hispanic), Ages 40-64
American Lung Association: https://www.lung.org/research/trends-in-lung-disease/tobacco-trends-brief/data-tables/ad-cig-smoke-rate-sex-race-age
- Percent of Adults Who Smoked >24 Cigarettes Daily by Sex, Race, Hispanic Origin, Age and Education
Dados <- readxl::read_excel("Lung_Sex.xlsx", sheet="LungMenSmoke24")
Dados <- timeSeries::timeSeries(data=Dados[,2:3],
charvec=timeDate::timeSequence(from="1997",
to="2018",
by="year"),
title="Age-adjusted incidence of Lung cancer
in USA")
Dados <- cbind(Dados,diff(log(Dados$Lung)))
Dados <- cbind(Dados,diff(log(Dados$Smoke)))
colnames(Dados) <- c("Lung", "Smoke", "LungDiffLn", "SmokeDiffLn")
print(Dados)GMT
Lung Smoke LungDiffLn SmokeDiffLn
1997-08-08 24.4 21.9 0.036221263 0.070513784
1998-08-08 25.3 23.5 -0.086652117 -0.151317817
1999-08-08 23.2 20.2 -0.076099344 0.024451096
2000-08-08 21.5 20.7 0.000000000 -0.106972120
2001-08-08 21.5 18.6 -0.009345862 -0.021739987
2002-08-08 21.3 18.2 -0.130183549 -0.122602322
2003-08-08 18.7 16.1 -0.119120828 0.024541109
2004-08-08 16.6 16.5 0.000000000 -0.006079046
2005-08-08 16.6 16.4 0.011976191 -0.062913825
2006-08-08 16.8 15.4 -0.018018506 -0.161755279
2007-08-08 16.5 13.1 0.018018506 -0.071176278
2008-08-08 16.8 12.2 0.040821995 -0.041847110
2009-08-08 17.5 11.7 -0.040821995 0.074107972
2010-08-08 16.8 12.6 0.133531393 0.068992871
2011-08-08 19.2 13.5 -0.064538521 -0.330563800
2012-08-08 18.0 9.7 0.011049836 0.060018010
2013-08-08 18.2 10.3 -0.068208250 0.038099846
2014-08-08 17.0 10.7 0.079137321 0.063369614
2015-08-08 18.4 11.4 -0.115069330 -0.082238098
2016-08-08 16.4 10.5 -0.024692613 -0.143100844
2017-08-08 16.0 9.1 -0.098440073 0.063851472
2018-08-08 14.5 9.7 0.036221263 0.070513784
plot(Dados$LungDiffLn, type="l", ylab="", ylim=c(-.3,.2))
lines(Dados$SmokeDiffLn, type="l", lty=2, lwd=2)
legend("topleft", legend=c("LungDiffLn","SunspotDiffLn"), lty=c(1,2), bty="n")
Pearson's product-moment correlation
data: Dados$LungDiffLn and Dados$SmokeDiffLn
t = 1.5662, df = 20, p-value = 0.133
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1058240 0.6601455
sample estimates:
cor
0.3305316 $index
integer(0)
$replacements
numeric(0)$index
integer(0)
$replacements
numeric(0)fitauto <- forecast::auto.arima(Dados$Lung,
stationary=TRUE,
ic="bic",
lambda=0, # Box-Cox transformation
xreg=Dados$Smoke)
summary(fitauto)Series: Dados$Lung
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 0
Coefficients:
ar1 intercept xreg
0.6609 2.5100 0.0277
s.e. 0.1687 0.1004 0.0064
sigma^2 = 0.004224: log likelihood = 30.25
AIC=-52.49 AICc=-50.14 BIC=-48.13
Training set error measures:
ME RMSE MAE MPE MAPE MASE
Training set -0.01944662 1.099527 0.9378088 -0.5107523 5.135195 0.8911306
ACF1
Training set 0.02769444
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 0.6608809 0.1686580 3.9185 8.911e-05 ***
intercept 2.5100412 0.1004481 24.9884 < 2.2e-16 ***
xreg 0.0276619 0.0063845 4.3327 1.473e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: fitauto$residuals
X-squared = 0.01887, df = 1, p-value = 0.8907Métodos de impactos de evento em série temporal
- Análise de quebra estrutural por SARIMA
- Análise de intervenção por SARIMAX
- Intervenção linear degrau permanente (step), puslo (impulso) e rampa (efeito tansiente)
- Box, GEP, Jenkins, GM, Reinsel, GC, and Ljung, GM (2015) Time series analysis: Forecasting and control [in R]. 5th ed. NJ: Wiley.
- Análise de função resposta ao impulso
- Análise de intervenção por GzLM em série temporal de contagem:
tscount::tsglmetscount::interv_test - Análise CAPM (beta)
- Valadkhani, A (2025) Inflation-driven instability in US sectoral betas. J Asset Manag (IF = 2.5). https://doi.org/10.1057/s41260-025-00413-3
- Análise de retorno anornal acumulado (CAR)
- Análise de volatilidade condicional por GARCH
- Análise de índice de Sharpe econométrico
Índices de bolsa de valores do setor de saúde (health care)
- Valadkhani, A (2025) Inflation-driven instability in US sectoral
betas. J Asset Manag (IF = 2.5). https://doi.org/10.1057/s41260-025-00413-3
- Resumo: Este estudo especifica e estima os coeficientes beta mensais como função da inflação anualizada para todos os fundos setoriais negociados em bolsa (ETF) nos EUA entre janeiro de 1999 e fevereiro de 2024. A consistência nas definições e a ausência de sobreposição de componentes são garantidas por meio do uso dos ETF setoriais da SPDR, que possuem a série histórica de retornos mais longa. Os resultados indicam que, quando a inflação sai da faixa de 1–5%, especialmente acima de 4–5%, os betas funcionais de todos os setores passam a apresentar comportamento inconsistente e tornam-se extremamente não confiáveis. Esse comportamento evidencia a vulnerabilidade de todos os ETFs setoriais à alta inflação. Estratégias de investimento devem considerar que a meta de inflação de 2% do Fed é crucial para garantir a confiabilidade do beta como medida de risco sistemático entre os setores. Esse achado ressalta, recentemente, a necessidade de o Fed controlar a inflação antes de reduzir as taxas de juros, visando à estabilidade dos mercados financeiros.
ETF (Exchange-Traded Fund) é um fundo de investimento negociado em bolsa, como se fosse uma ação. Ele replica o desempenho de um índice, setor, commodity ou outro ativo, permitindo ao investidor obter diversificação automática.
Exemplos:
- SPY: replica o S&P 500
- XLV: setor de saúde dos EUA
- QQQ: Nasdaq-100
- EWW: ações do México
Características principais:
- Negociação em tempo real (como ações)
- Taxas de administração geralmente baixas
- Liquidez alta (dependendo do ativo)
- Pode ser usado para hedge, alocação estratégica ou exposição tática
No contexto do CAPM, ETF são úteis para medir o beta setorial e analisar o risco sistemático relativo de diferentes setores.
XLV é um ETF da Health Care Select Sector SPDR Fund que replica o desempenho do setor de saúde do S&P 500. Inclui empresas farmacêuticas, de biotecnologia, equipamentos médicos e operadoras de planos de saúde. É considerado um setor defensivo, com estimativa pontual de beta historicamente abaixo de 1 (0.76 (2018) caindo para 0.71 (2024)) (ver Fig. 2).
Outros indices do setor de Health Care dos EUA:
- Nasdaq Health Care Index
- Dow Jones U.S. Health Care Index
- MSCI U.S. IMI Health Care 25/50 Index
- NYSE Arca Biotechnology Index
Valadkhani, 2025, Table 2
Valadkhani, 2025, Fig. 2
Análise de intervenção de COVID-19 em índice de bolsa de valores: SARIMAX
SARIMAX
SARIMAX (Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors) é uma extensão do modelo ARIMA que incorpora:
Sazonalidade (SARIMA) Modela padrões que se repetem ciclicamente (ex: mensal, semanal).
Variáveis Exógenas (X) Permite incluir covariáveis externas que afetam a série temporal (ex: políticas públicas, pandemias, choques econômicos).
Formalmente:
\[ Y_t = \text{SARIMA}(p,d,q)(P,D,Q)_s + \boldsymbol{\beta}^{\prime} \boldsymbol{X}_t + \varepsilon_t \]
Onde:
- \(p, d, q\): ordens do ARIMA
- \(P, D, Q\): ordens sazonais
- \(s\): periodicidade (ex: \(s = 12\) para mensal)
- \(\boldsymbol{X}_t\): vetor de variáveis exógenas
- \(\boldsymbol{\beta}\): coeficientes de regressão
- \(\varepsilon_t \sim \text{Ruído branco}\)
Em R, pode ser ajustado via:
forecast::Arima(y,
order = c(p,d,q),
seasonal = list(order = c(P,D,Q), period = s),
xreg = X)
Uso típico em saúde ou economia: avaliar o impacto de intervenções (ex: pandemia, vacinação) com controle de sazonalidade.
Modelo de intervenção de degrau
Modelo de intervenção tipo degrau (step) serve para medir o efeito de um evento que causa uma mudança permanente no nível da série temporal a partir de uma data específica.
Exemplo: impacto da pandemia de COVID-19 no número de consultas médicas.
Funciona assim:
- Antes da data do evento, o modelo estima a série como ela vinha normalmente.
- A partir da data do evento, ele adiciona um efeito constante (positivo ou negativo), representando a nova realidade da série.
- Esse efeito é somado ao modelo ARIMA que já explica a série antes do evento.
A variável “degrau” é binária:
- 0 antes do evento
- 1 durante o evento
- 0 após o evento
O modelo estima quanto a média da série mudou após o evento. É útil quando há mudanças estruturais permanentes, como políticas públicas, crises ou intervenções de saúde.
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Baixar dados do IBOVESPA
quantmod::getSymbols("^BVSP", from = "2018-01-01", to = "2024-12-31", src = "yahoo")[1] "BVSP"ibov <- zoo::na.locf(quantmod::Cl(BVSP))
ibov_lr <- diff(log(ibov))
# Datas dos eventos
inicio_pandemia <- as.Date("2020-03-16")
inicio_vacina <- as.Date("2021-01-18")
fim_pandemia <- as.Date("2023-05-05")
# Indexação da série
dias_index <- zoo::index(ibov_lr)
n <- length(dias_index)
# Vetores de intervenção (step functions)
step_semvacina <- ifelse(dias_index >= inicio_pandemia & dias_index < inicio_vacina, 1, 0)
step_comvacina <- ifelse(dias_index >= inicio_vacina & dias_index <= fim_pandemia, 1, 0)
# Matriz de regressoras exógenas
X_step <- cbind(step_semvacina, step_comvacina)
# Data frame para gráfico
df_step <- data.frame(
Data = dias_index,
Step_sem_vacina = step_semvacina,
Step_com_vacina = step_comvacina
)
# Gráfico
g <- ggplot2::ggplot(df_step, ggplot2::aes(x = Data)) +
ggplot2::geom_step(ggplot2::aes(y = Step_sem_vacina, color = "Sem vacina")) +
ggplot2::geom_step(ggplot2::aes(y = Step_com_vacina, color = "Com vacina")) +
ggplot2::scale_color_manual(values = c("Sem vacina" = "red", "Com vacina" = "blue")) +
ggplot2::labs(title = "Funções de Intervenção (Step) na Pandemia de COVID-19",
y = "Valor da Intervenção", x = "Data", color = "Fase") +
ggplot2::theme_minimal(base_size = 14)
plot(g)
# Ajustar modelo ARIMA com step
modelo_step <- forecast::auto.arima(ibov_lr,
stationary = TRUE,
ic = "bic",
lambda = 1,
xreg = X_step)
# Diagnóstico do modelo
print(summary(modelo_step))Series: ibov_lr
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept step_semvacina step_comvacina
-0.1437 -0.9998 0.0018 -4e-04
s.e. 0.0237 0.0004 0.0010 7e-04
sigma^2 = 0.0002341: log likelihood = 4797.79
AIC=-9585.58 AICc=-9585.55 BIC=-9558.28
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 6.030963e-07 0.01528233 0.01043279 Inf Inf 0.6850147 0.006013344
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.14374940 0.02374754 -6.0532 1.42e-09 ***
intercept -0.99982814 0.00043271 -2310.6339 < 2.2e-16 ***
step_semvacina 0.00180560 0.00102335 1.7644 0.07766 .
step_comvacina -0.00042039 0.00070678 -0.5948 0.55198
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: modelo_step$residuals
X-squared = 0.06281, df = 1, p-value = 0.8021# Ajustar modelo ARIMA(1,0,0) com regressoras
modelo_step <- forecast::Arima(ibov_lr,
order = c(1, 0, 0),
xreg = X_step)
# Resumo e diagnóstico
summary(modelo_step)Series: ibov_lr
Regression with ARIMA(1,0,0) errors
Coefficients:
ar1 intercept step_semvacina step_comvacina
-0.1437 2e-04 0.0018 -4e-04
s.e. 0.0237 4e-04 0.0010 7e-04
sigma^2 = 0.0002341: log likelihood = 4797.79
AIC=-9585.58 AICc=-9585.55 BIC=-9558.28
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 7.292659e-08 0.01528233 0.01043278 Inf Inf 0.6852321 0.006010624
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -0.14374663 0.02374754 -6.0531 1.421e-09 ***
intercept 0.00017264 0.00043271 0.3990 0.68992
step_semvacina 0.00180518 0.00102335 1.7640 0.07774 .
step_comvacina -0.00042119 0.00070678 -0.5959 0.55122
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: modelo_step$residuals
X-squared = 0.062754, df = 1, p-value = 0.8022# Gráfico final
plot(zoo::index(ibov), ibov, type = "l", col = "blue", lwd = 2,
main = "Impacto da Pandemia de COVID-19 no IBOV\nModelo Step",
xlab = "Tempo", ylab = "Fechamento")
abline(v = inicio_pandemia, col = "red", lty = 2, lwd = 2)
abline(v = inicio_vacina, col = "darkgreen", lty = 2, lwd = 2)
abline(v = fim_pandemia, col = "purple", lty = 2, lwd = 2)
legend("bottomright",
legend = c("IBOV", "Início pandemia", "Início vacinação", "Fim pandemia"),
col = c("blue", "red", "darkgreen", "purple"),
lty = c(1, 2, 2, 2),
lwd = c(2, 2, 2, 2),
bty = "n")
Modelo de intervenção linear de rampa
O modelo de intervenção por rampa linear crescente e decrescente é usado para representar mudanças graduais em séries temporais após eventos que não provocam efeitos imediatos, mas sim acumulativos ou reversíveis ao longo do tempo.
No contexto da pandemia de COVID-19, para um índice como o IBOV ou XLV, definimos duas fases de intervenção:
Rampa linear crescente (sem vacina): Começa no início da pandemia (ex.: 2020-03-16) e aumenta linearmente a cada dia até o início da vacinação (ex.: 2021-01-17). Captura o acúmulo de impacto negativo da pandemia sobre o mercado.
\[ R_1(t) = \begin{cases} t - t_0 + 1, & \text{se } t_0 \le t < t_1 \\ 0, & \text{caso contrário} \end{cases} \]
Rampa linear decrescente (com vacina): Começa no início da vacinação (ex.: 2021-01-18) com valor igual ao final da rampa crescente, e decresce linearmente até 0 ao final da pandemia (ex.: 2023-05-04). Representa a redução gradual dos efeitos da pandemia.
\[ R_2(t) = \begin{cases} L - (t - t_1), & \text{se } t_1 \le t \le t_2 \\ 0, & \text{caso contrário} \end{cases} \]
Onde \(L = t_1 - t_0\) é o comprimento da rampa crescente.
Essas duas variáveis de intervenção são incluídas como regressores exógenos (xreg) em modelos ARIMA ou regressão com erros autocorrelacionados, para estimar e testar o impacto gradual e reversível da pandemia sobre a série.
Esse tipo de modelagem é mais realista do que uma variável binária, pois reflete a dinâmica não instantânea dos efeitos econômicos.
IBOVESPA
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Baixar dados do IBOVESPA até final de 2024
quantmod::getSymbols("^BVSP", from = "2018-01-01", to = "2024-12-31", src = "yahoo")[1] "BVSP"ibov <- zoo::na.locf(quantmod::Cl(BVSP)) # Preço de fechamento ajustado
ibov_lr <- diff(log(ibov))
# Datas das fases
inicio_pandemia <- as.Date("2020-03-16")
inicio_vacina <- as.Date("2021-01-18")
fim_pandemia <- as.Date("2023-05-05") # OMS declara fim
# Gráfico final com três marcos temporais (sem linha de ajuste)
plot(zoo::index(ibov), ibov, type = "l", col = "blue", lwd = 2,
main = "Impacto da Pandemia de COVID-19 no IBOV",
xlab = "Tempo", ylab = "Fechamento")
abline(v = inicio_pandemia, col = "red", lty = 2, lwd = 2)
abline(v = inicio_vacina, col = "darkgreen", lty = 2, lwd = 2)
abline(v = fim_pandemia, col = "purple", lty = 2, lwd = 2)
legend("topleft",
legend = c("IBOV", "Início pandemia", "Início vacinação", "Fim pandemia"),
col = c("blue", "red", "darkgreen", "purple"),
lty = c(1, 2, 2, 2),
lwd = c(2, 2, 2, 2),
bty = "n")
# Converter para vetor e objeto ts
ibov_vec <- as.numeric(ibov)
ibov_ts <- stats::ts(ibov_vec, start = c(2018, 1), frequency = 252)
Dados <- ibov_ts
cpm <- changepoint::cpt.mean(Dados)
print(summary(cpm))Created Using changepoint version 2.3
Changepoint type : Change in mean
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 22.38147
Minimum Segment Length : 1
Maximum no. of cpts : Inf
Number of changepoints: 1732
NULL
Created Using changepoint version 2.3
Changepoint type : Change in variance
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 22.38147
Minimum Segment Length : 2
Maximum no. of cpts : Inf
Changepoint Locations : 246 361 447 477 539 598 810 918 1009 1207 1345 1456
NULL
Optimal 5-segment partition:
Call:
breakpoints.formula(formula = Dados ~ 1)
Breakpoints at observation number:
260 712 972 1454
Corresponding to breakdates:
2019(8) 2020(208) 2021(216) 2023(194)
Optimal (m+1)-segment partition:
Call:
breakpoints.formula(formula = Dados ~ 1)
Breakpoints at observation number:
m = 1 712
m = 2 707 1456
m = 3 260 712 1456
m = 4 260 712 972 1454
m = 5 260 694 954 1214 1474
Corresponding to breakdates:
m = 1 2020(208)
m = 2 2020(203) 2023(196)
m = 3 2019(8) 2020(208) 2023(196)
m = 4 2019(8) 2020(208) 2021(216) 2023(194)
m = 5 2019(8) 2020(190) 2021(198) 2022(206) 2023(214)
Fit:
m 0 1 2 3 4 5
RSS 4.384e+11 1.909e+11 1.427e+11 9.713e+10 8.990e+10 9.516e+10
BIC 3.857e+04 3.714e+04 3.665e+04 3.600e+04 3.588e+04 3.599e+04
# Indexação da série
dias_index <- zoo::index(ibov)
n <- length(dias_index)
# Vetores das rampas
interv_semvacina <- rep(0, n)
interv_comvacina <- rep(0, n)
# Índices das fases
idx_semvacina <- which(dias_index >= inicio_pandemia & dias_index < inicio_vacina)
idx_comvacina <- which(dias_index >= inicio_vacina & dias_index <= fim_pandemia)
# Rampa crescente: pandemia sem vacina
interv_semvacina[idx_semvacina] <- seq_len(length(idx_semvacina))
# Rampa decrescente: pandemia com vacina até 0 no fim
valor_inicial <- max(interv_semvacina) # agora pega corretamente o nível final
interv_comvacina[idx_comvacina] <- seq(from = valor_inicial,
to = 0,
length.out = length(idx_comvacina))
# Data frame para plotagem
df_rampas <- data.frame(
Data = dias_index,
Rampa_sem_vacina = interv_semvacina,
Rampa_com_vacina = interv_comvacina
)
# Plotar as duas rampas
g <- ggplot2::ggplot(df_rampas, ggplot2::aes(x = Data)) +
ggplot2::geom_line(ggplot2::aes(y = Rampa_sem_vacina, color = "Sem vacina")) +
ggplot2::geom_line(ggplot2::aes(y = Rampa_com_vacina, color = "Com vacina")) +
ggplot2::scale_color_manual(values = c("Sem vacina" = "red", "Com vacina" = "blue")) +
ggplot2::labs(title = "Rampas de Intervenção na Pandemia",
y = "Nível da intervenção", x = "Data", color = "Fase") +
ggplot2::theme_minimal(base_size = 14)
plot(g)
# Combinar como matriz
X <- cbind(sem_vacina = interv_semvacina,
com_vacina = interv_comvacina)
# Ajustar modelo ARIMA com duas rampas
modelo_duplo <- forecast::auto.arima(ibov_lr,
stationary = TRUE,
ic = "bic",
lambda = 1,
xreg = X)
# Diagnóstico do modelo
print(summary(modelo_duplo))Series: ibov_lr
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept sem_vacina com_vacina
-0.1429 -0.9998 0 0
s.e. 0.0237 0.0004 0 0
sigma^2 = 0.0002341: log likelihood = 4797.73
AIC=-9585.46 AICc=-9585.42 BIC=-9558.16
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 6.175275e-07 0.01528288 0.01043566 Inf Inf 0.685203 0.006024488
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -1.4285e-01 2.3745e-02 -6.0161 1.787e-09 ***
intercept -9.9983e-01 3.9477e-04 -2532.7041 < 2.2e-16 ***
sem_vacina 1.4887e-05 2.5359e-05 0.5871 0.5572
com_vacina -3.2402e-06 2.4606e-05 -0.1317 0.8952
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: modelo_duplo$residuals
X-squared = 0.063043, df = 1, p-value = 0.8017# Ajustar modelo ARIMA com duas rampas
modelo_duplo <- forecast::Arima(ibov_lr,
order=c(1,0,0),
lambda = 1,
xreg = X)
# Diagnóstico do modelo
print(summary(modelo_duplo))Series: ibov_lr
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept sem_vacina com_vacina
-0.1429 -0.9998 0 0
s.e. 0.0237 0.0004 0 0
sigma^2 = 0.0002341: log likelihood = 4797.73
AIC=-9585.46 AICc=-9585.42 BIC=-9558.16
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 8.595558e-08 0.01528288 0.0104357 Inf Inf 0.685424 0.006090677
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -1.4292e-01 2.3745e-02 -6.0187 1.759e-09 ***
intercept -9.9982e-01 3.9474e-04 -2532.8420 < 2.2e-16 ***
sem_vacina 1.4874e-05 2.5359e-05 0.5866 0.5575
com_vacina -3.2425e-06 2.4606e-05 -0.1318 0.8952
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: modelo_duplo$residuals
X-squared = 0.064436, df = 1, p-value = 0.7996ETF XLV
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Baixar dados do XLV até final de 2024
quantmod::getSymbols("XLV", from = "2018-01-01", to = "2024-12-31", src = "yahoo")[1] "XLV"xlv <- zoo::na.locf(quantmod::Cl(XLV)) # Preço de fechamento ajustado
xlv_lr <- diff(log(xlv))
# Datas das fases
inicio_pandemia <- as.Date("2020-03-16")
inicio_vacina <- as.Date("2021-01-18")
fim_pandemia <- as.Date("2023-05-05") # OMS declara fim
# Gráfico final com três marcos temporais (sem linha de ajuste)
plot(zoo::index(xlv), xlv, type = "l", col = "blue", lwd = 2,
main = "Impacto da Pandemia de COVID-19 no XLV",
xlab = "Tempo", ylab = "Fechamento")
abline(v = inicio_pandemia, col = "red", lty = 2, lwd = 2)
abline(v = inicio_vacina, col = "darkgreen", lty = 2, lwd = 2)
abline(v = fim_pandemia, col = "purple", lty = 2, lwd = 2)
legend("topleft",
legend = c("xlv", "Início pandemia", "Início vacinação", "Fim pandemia"),
col = c("blue", "red", "darkgreen", "purple"),
lty = c(1, 2, 2, 2),
lwd = c(2, 2, 2, 2),
bty = "n")
# Converter para vetor e objeto ts
xlv_vec <- as.numeric(xlv)
xlv_ts <- stats::ts(xlv_vec, start = c(2018, 1), frequency = 252)
Dados <- xlv_ts
cpm <- changepoint::cpt.mean(Dados)
print(summary(cpm))Created Using changepoint version 2.3
Changepoint type : Change in mean
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 22.41921
Minimum Segment Length : 1
Maximum no. of cpts : Inf
Number of changepoints: 92
NULL
Created Using changepoint version 2.3
Changepoint type : Change in variance
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 22.41921
Minimum Segment Length : 2
Maximum no. of cpts : Inf
Changepoint Locations : 484 638 755 826 880 1511
NULL
Optimal 5-segment partition:
Call:
breakpoints.formula(formula = Dados ~ 1)
Breakpoints at observation number:
271 575 839 1496
Corresponding to breakdates:
2019(19) 2020(71) 2021(83) 2023(236)
Optimal (m+1)-segment partition:
Call:
breakpoints.formula(formula = Dados ~ 1)
Breakpoints at observation number:
m = 1 758
m = 2 735 1496
m = 3 472 822 1496
m = 4 271 575 839 1496
m = 5 271 574 838 1213 1496
Corresponding to breakdates:
m = 1 2021(2)
m = 2 2020(231) 2023(236)
m = 3 2019(220) 2021(66) 2023(236)
m = 4 2019(19) 2020(71) 2021(83) 2023(236)
m = 5 2019(19) 2020(70) 2021(82) 2022(205) 2023(236)
Fit:
m 0 1 2 3 4 5
RSS 818937 138428 85554 49456 44398 44036
BIC 15821 12707 11875 10925 10750 10751
# Indexação da série
dias_index <- zoo::index(xlv)
n <- length(dias_index)
# Vetores das rampas
interv_semvacina <- rep(0, n)
interv_comvacina <- rep(0, n)
# Índices das fases
idx_semvacina <- which(dias_index >= inicio_pandemia & dias_index < inicio_vacina)
idx_comvacina <- which(dias_index >= inicio_vacina & dias_index <= fim_pandemia)
# Rampa crescente: pandemia sem vacina
interv_semvacina[idx_semvacina] <- seq_len(length(idx_semvacina))
# Rampa decrescente: pandemia com vacina até 0 no fim
valor_inicial <- max(interv_semvacina) # agora pega corretamente o nível final
interv_comvacina[idx_comvacina] <- seq(from = valor_inicial,
to = 0,
length.out = length(idx_comvacina))
# Data frame para plotagem
df_rampas <- data.frame(
Data = dias_index,
Rampa_sem_vacina = interv_semvacina,
Rampa_com_vacina = interv_comvacina
)
# Plotar as duas rampas
g <- ggplot2::ggplot(df_rampas, ggplot2::aes(x = Data)) +
ggplot2::geom_line(ggplot2::aes(y = Rampa_sem_vacina, color = "Sem vacina")) +
ggplot2::geom_line(ggplot2::aes(y = Rampa_com_vacina, color = "Com vacina")) +
ggplot2::scale_color_manual(values = c("Sem vacina" = "red", "Com vacina" = "blue")) +
ggplot2::labs(title = "Rampas de Intervenção na Pandemia",
y = "Nível da intervenção", x = "Data", color = "Fase") +
ggplot2::theme_minimal(base_size = 14)
plot(g)
# Combinar como matriz
X <- cbind(sem_vacina = interv_semvacina,
com_vacina = interv_comvacina)
# Ajustar modelo ARIMA com duas rampas
modelo_duplo <- forecast::auto.arima(xlv_lr,
stationary = TRUE,
ic = "bic",
lambda = 1,
xreg = X)
# Diagnóstico do modelo
print(summary(modelo_duplo))Series: xlv_lr
Regression with ARIMA(2,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 ar2 intercept sem_vacina com_vacina
-0.0958 0.0930 -0.9999 0 0
s.e. 0.0237 0.0237 0.0003 0 0
sigma^2 = 0.0001214: log likelihood = 5436.58
AIC=-10861.17 AICc=-10861.12 BIC=-10828.33
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 8.200015e-07 0.01100208 0.007600622 Inf Inf 0.7040351 0.001326591
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -9.5805e-02 2.3749e-02 -4.0341 5.482e-05 ***
ar2 9.3019e-02 2.3742e-02 3.9180 8.930e-05 ***
intercept -9.9986e-01 3.2230e-04 -3102.2632 < 2.2e-16 ***
sem_vacina 7.1447e-06 2.4713e-05 0.2891 0.7725
com_vacina 1.4053e-06 2.4236e-05 0.0580 0.9538
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: modelo_duplo$residuals
X-squared = 0.0030956, df = 1, p-value = 0.9556# Ajustar modelo ARIMA com duas rampas
modelo_duplo <- forecast::Arima(xlv_lr,
order=c(1,0,0),
lambda = 1,
xreg = X)
# Diagnóstico do modelo
print(summary(modelo_duplo))Series: xlv_lr
Regression with ARIMA(1,0,0) errors
Box Cox transformation: lambda= 1
Coefficients:
ar1 intercept sem_vacina com_vacina
-0.1057 -0.9999 0 0
s.e. 0.0237 0.0003 0 0
sigma^2 = 0.0001224: log likelihood = 5428.94
AIC=-10847.88 AICc=-10847.85 BIC=-10820.52
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 6.840508e-07 0.01105004 0.007600588 Inf Inf 0.7041463 0.009815923
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
ar1 -1.0565e-01 2.3719e-02 -4.4542 8.42e-06 ***
intercept -9.9986e-01 2.9380e-04 -3403.1839 < 2.2e-16 ***
sem_vacina 7.1502e-06 2.4568e-05 0.2910 0.7710
com_vacina 1.3872e-06 2.4169e-05 0.0574 0.9542
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Box-Pierce test
data: modelo_duplo$residuals
X-squared = 0.16948, df = 1, p-value = 0.6806Análise com CAPM do efeito de COVID-19 em ETF XLV
Risco sistemático
risco sistemático (ou não diversificável) é o risco inerente ao mercado como um todo, que afeta todos os ativos simultaneamente, independentemente das características individuais.
Exemplos típicos:
- Crises econômicas
- Pandemias (ex.: COVID-19)
- Choques geopolíticos
- Políticas monetárias
No CAPM, o risco sistemático é medido pelo beta (\(\beta\)):
- \(\beta = 1\): risco igual ao mercado
- \(\beta < 1\): risco menor que o mercado (ativo defensivo)
- \(\beta > 1\): risco maior que o mercado (ativo agressivo)
risco sistemático não pode ser eliminado por diversificação. Já o risco específico (ou idiossincrático) pode ser eliminado ao combinar ativos distintos em carteira.
CAPM
O modelo CAPM (Capital Asset Pricing Model) expressa o retorno esperado de um ativo em função de seu risco sistemático:
\[ \mathbb{E}(R_i) = R_f + \beta_i \, (\mathbb{E}(R_m) - R_f) \]
Onde:
- \(\mathbb{E}(R_i)\): retorno esperado do ativo \(i\)
- \(R_f\): taxa livre de risco
- \(\mathbb{E}(R_m)\): retorno esperado do mercado
- \(\beta_i\): sensibilidade do ativo ao mercado (risco sistemático)
Na prática, estima-se via regressão linear:
\[ (R_i - R_f) = \alpha + \beta (R_m - R_f) + \varepsilon \]
Ou seja, retorno em excesso do ativo sobre a taxa livre de risco é explicado pelo retorno em excesso do mercado.
\(\beta\) indica o quanto o ativo acompanha o mercado.
\(\alpha\) positivo sugere retorno acima do esperado (\(\alpha\) anômalo).
Beta do CAPM
O beta do CAPM é uma medida de risco sistemático — ou seja, risco não diversificável — que indica a sensibilidade de um ativo (ou setor) às variações do mercado.
Para o setor de saúde (ex.: ETF XLV):
- Se \(\beta < 1\), o setor tende a reagir menos do que o mercado a choques sistêmicos (setor defensivo).
- Se \(\beta > 1\), é mais volátil que o mercado (não é o caso do XLV historicamente).
- Exemplo: \(\beta = 0.65\) implica que uma alta de 1% no mercado gera, em média, alta de 0.65% no XLV.
Na pandemia, a queda progressiva dos betas do XLV sugere que o setor ganhou caráter mais defensivo, absorvendo menor risco sistemático, possivelmente devido à sua relevância sanitária e estabilidade de receitas.
Dia útil
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Baixar dados de XLV, S&P500 (GSPC) e taxa livre de risco (IRX = 13-week T-Bill)
quantmod::getSymbols(c("XLV", "^GSPC", "^IRX"), from = "2018-01-01", to = "2023-12-31")[1] "XLV" "GSPC" "IRX" # Calcular retornos logarítmicos diários
r_XLV <- diff(log(quantmod::Cl(XLV)))
r_GSPC <- diff(log(quantmod::Cl(GSPC)))
rf <- (quantmod::Cl(IRX) / 100 + 1)^(1/252) - 1 # taxa livre de risco ao dia util
# Alinhar e calcular retornos em excesso
dados <- na.omit(merge(r_XLV, r_GSPC, rf))
colnames(dados) <- c("r_XLV", "r_GSPC", "r_f")
dados$excess_XLV <- dados$r_XLV - dados$r_f
dados$excess_GSPC <- dados$r_GSPC - dados$r_f
# Definir os períodos
periodos <- list(
pre_pandemia = c("2018-01-01", "2020-03-15"),
sem_vacina = c("2020-03-16", "2021-01-17"),
com_vacina = c("2021-01-18", "2023-05-04"),
pos_pandemia = c("2023-05-05", "2024-12-31")
)
alfa <- 0.05
# Estimar modelos e ICs com Newey-West
resultados_nw <- lapply(names(periodos), function(nome) {
intervalo <- periodos[[nome]]
dados_fase <- window(dados, start = as.Date(intervalo[1]), end = as.Date(intervalo[2]))
modelo <- lm(excess_XLV ~ excess_GSPC, data = dados_fase)
# Estimativa robusta do erro padrão com Newey-West
nw_se <- sqrt(diag(sandwich::NeweyWest(modelo,
lag = 5,
prewhite = FALSE)))
est <- coef(modelo)["excess_GSPC"]
tcrit <- qt(1-alfa/4, df = modelo$df.residual)
ic_inf <- est - tcrit * nw_se["excess_GSPC"]
ic_sup <- est + tcrit * nw_se["excess_GSPC"]
data.frame(Periodo = nome,
Beta = est,
IC_inferior = ic_inf,
IC_superior = ic_sup)
})
# Tabela final
df_nw <- do.call(rbind, resultados_nw)
print(df_nw, digits = 3) Periodo Beta IC_inferior IC_superior
excess_GSPC pre_pandemia 0.853 0.777 0.930
excess_GSPC1 sem_vacina 0.821 0.756 0.886
excess_GSPC2 com_vacina 0.633 0.586 0.680
excess_GSPC3 pos_pandemia 0.565 0.441 0.689# Garantir ordem dos fatores
df_nw$Periodo <- factor(df_nw$Periodo,
levels = c("pre_pandemia", "sem_vacina",
"com_vacina", "pos_pandemia"))
# Gráfico com sintaxe explícita
g <- ggplot2::ggplot(df_nw, ggplot2::aes(x = Periodo, y = Beta)) +
ggplot2::geom_point(size = 3, color = "blue") +
ggplot2::geom_errorbar(ggplot2::aes(ymin = IC_inferior, ymax = IC_superior),
width = 0.2, color = "gray40") +
ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
ggplot2::labs(title = "Beta do CAPM para XLV\nIC95% Bonf robusto (Newey-West)",
x = "Período", y = "Beta") +
ggplot2::theme_minimal(base_size = 14)
plot(g)
A interpretação dos betas do CAPM estimados para o ETF XLV (setor de saúde dos EUA) em diferentes fases da pandemia de COVID-19 é a seguinte:
Pré-pandemia (β = 0,853; IC95%: 0,777 a 0,930): O setor de saúde era moderadamente sensível ao mercado. Uma variação de 1% no S&P500 resultava, em média, em uma variação de 0,85% no XLV. O risco sistemático era significativamente diferente de 0 e menor que 1 (defensivo).
Pandemia sem vacina (β = 0,821; IC95%: 0,756 a 0,886): Houve leve queda no beta, sugerindo estabilidade do perfil defensivo mesmo sob incerteza extrema. O intervalo de confiança se sobrepõe ao pré-pandemia ⇒ sem evidência forte de mudança estatística.
Pandemia com vacina (β = 0,633; IC95%: 0,586 a 0,680): Redução clara do beta. O setor ficou ainda menos sensível ao mercado. A vacinação pode ter aumentado a previsibilidade dos lucros e reduzido o risco sistemático. IC não se sobrepõe ao período anterior ⇒ evidência de mudança.
Pós-pandemia (β = 0,565; IC95%: 0,441 a 0,689): Beta continua baixo. Setor manteve comportamento defensivo. Intervalo de confiança se sobrepõe ao da fase com vacina ⇒ sem evidência forte de nova mudança.
Resumo: O setor de saúde (XLV) manteve-se defensivo em todas as fases, com redução significativa do risco sistemático após o início da vacinação. O beta < 1 reforça o caráter resiliente do setor frente às flutuações do mercado.
Semanal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Baixar dados de preços
quantmod::getSymbols(c("XLV", "^GSPC", "^IRX"), from = "2018-01-01", to = "2023-12-31")[1] "XLV" "GSPC" "IRX" # Retornos semanais logarítmicos
r_XLV <- quantmod::weeklyReturn(XLV, type = "log")
r_GSPC <- quantmod::weeklyReturn(GSPC, type = "log")
# Taxa livre de risco semanal: IRX/100/52 (conversão de % anual para semanal)
rf_semanal <- zoo::na.locf(quantmod::Cl(IRX)) / 100 / 52
rf_semanal <- xts::apply.weekly(rf_semanal, mean)NOTE: `apply.weekly(..., FUN = mean)` operates by column, unlike other math
functions (e.g. median, sum, var, sd). Please use `FUN = colMeans` instead,
and use `FUN = function(x) mean(x)` to take the mean of all columns. Set
`options(xts.message.period.apply.mean = FALSE)` to suppress this message.# Dataset com retornos em excesso
dados <- na.omit(merge(r_XLV, r_GSPC, rf_semanal))
colnames(dados) <- c("r_XLV", "r_GSPC", "r_f")
dados$excess_XLV <- dados$r_XLV - dados$r_f
dados$excess_GSPC <- dados$r_GSPC - dados$r_f
# Períodos da pandemia
periodos <- list(
pre_pandemia = c("2018-01-01", "2020-03-15"),
sem_vacina = c("2020-03-16", "2021-01-17"),
com_vacina = c("2021-01-18", "2023-05-04"),
pos_pandemia = c("2023-05-05", "2023-12-31")
)
# Correção de Bonferroni para 4 comparações
alfa <- 0.05
# Estimar CAPM com broom::tidy() para cada fase
resultados <- lapply(names(periodos), function(nome) {
intervalo <- periodos[[nome]]
dados_fase <- window(dados, start = as.Date(intervalo[1]), end = as.Date(intervalo[2]))
modelo <- lm(excess_XLV ~ excess_GSPC, data = dados_fase)
resumo <- broom::tidy(modelo, conf.int = TRUE, conf.level = 1 - alfa/4)
beta <- resumo[resumo$term == "excess_GSPC", ]
data.frame(Periodo = nome,
Beta = beta$estimate,
IC_inferior = beta$conf.low,
IC_superior = beta$conf.high)
})
# Tabela final
df_resultados <- do.call(rbind, resultados)
df_resultados$Periodo <- factor(df_resultados$Periodo,
levels = c("pre_pandemia", "sem_vacina", "com_vacina", "pos_pandemia"))
print(df_resultados, digits = 3) Periodo Beta IC_inferior IC_superior
1 pre_pandemia 0.984 0.879 1.088
2 sem_vacina 0.811 0.617 1.006
3 com_vacina 0.639 0.508 0.770
4 pos_pandemia 0.663 0.347 0.979# Gráfico com ICs
g <- ggplot2::ggplot(df_resultados, ggplot2::aes(x = Periodo, y = Beta)) +
ggplot2::geom_point(size = 3, color = "blue") +
ggplot2::geom_errorbar(ggplot2::aes(ymin = IC_inferior, ymax = IC_superior),
width = 0.2, color = "darkgray") +
ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
ggplot2::labs(title = "Betas do CAPM para XLV por Fase da Pandemia\nRetorno logarítmico semanal",
x = "Período", y = "Beta") +
ggplot2::theme_minimal(base_size = 14)
plot(g)
Mensal
invisible(Sys.setlocale("LC_CTYPE", "pt_BR.UTF-8"))
invisible(Sys.setlocale("LC_ALL", "pt_BR.UTF-8"))
# Baixar dados
quantmod::getSymbols(c("XLV", "^GSPC", "^IRX"), from = "2018-01-01", to = "2024-12-31")[1] "XLV" "GSPC" "IRX" # Retornos logarítmicos mensais
r_XLV <- quantmod::monthlyReturn(XLV, type = "log")
r_GSPC <- quantmod::monthlyReturn(GSPC, type = "log")
# Taxa livre de risco convertida de anual para mensal
rf_mensal <- (zoo::na.locf(quantmod::Cl(IRX)) / 100 + 1)^(1/12) - 1
rf_mensal <- xts::apply.monthly(rf_mensal, mean)NOTE: `apply.monthly(..., FUN = mean)` operates by column, unlike other math
functions (e.g. median, sum, var, sd). Please use `FUN = colMeans` instead,
and use `FUN = function(x) mean(x)` to take the mean of all columns. Set
`options(xts.message.period.apply.mean = FALSE)` to suppress this message.# Combinar e calcular retornos em excesso
dados <- na.omit(merge(r_XLV, r_GSPC, rf_mensal))
colnames(dados) <- c("r_XLV", "r_GSPC", "r_f")
dados$excess_XLV <- dados$r_XLV - dados$r_f
dados$excess_GSPC <- dados$r_GSPC - dados$r_f
# Definir os períodos
periodos <- list(
pre_pandemia = c("2018-01-01", "2020-03-15"),
sem_vacina = c("2020-03-16", "2021-01-17"),
com_vacina = c("2021-01-18", "2023-05-04"),
pos_pandemia = c("2023-05-05", "2024-12-31")
)
alfa <- 0.05
# Estimar CAPM com broom::tidy()
resultados <- lapply(names(periodos), function(nome) {
intervalo <- periodos[[nome]]
dados_fase <- window(dados, start = as.Date(intervalo[1]), end = as.Date(intervalo[2]))
modelo <- lm(excess_XLV ~ excess_GSPC, data = dados_fase)
resumo <- broom::tidy(modelo, conf.int = TRUE, conf.level = 1 - alfa/4)
beta <- resumo[resumo$term == "excess_GSPC", ]
data.frame(Periodo = nome,
Beta = beta$estimate,
IC_inferior = beta$conf.low,
IC_superior = beta$conf.high)
})
# Combinar resultados
df_resultados <- do.call(rbind, resultados)
df_resultados$Periodo <- factor(df_resultados$Periodo,
levels = c("pre_pandemia", "sem_vacina", "com_vacina", "pos_pandemia"))
print(df_resultados, digits = 2) Periodo Beta IC_inferior IC_superior
1 pre_pandemia 0.85 0.52 1.17
2 sem_vacina 0.63 0.24 1.03
3 com_vacina 0.65 0.37 0.94
4 pos_pandemia 0.81 0.39 1.24# Gráfico
g <- ggplot2::ggplot(df_resultados, ggplot2::aes(x = Periodo, y = Beta)) +
ggplot2::geom_point(size = 3, color = "blue") +
ggplot2::geom_errorbar(ggplot2::aes(ymin = IC_inferior, ymax = IC_superior),
width = 0.2, color = "darkgray") +
ggplot2::geom_hline(yintercept = 1, linetype = "dashed", color = "red") +
ggplot2::labs(title = "Betas do CAPM por Fase da Pandemia (XLV)\nRetorno em excesso mensal",
x = "Período", y = "Beta") +
ggplot2::theme_minimal(base_size = 14)
plot(g)
Análise de intervenção por GzLM de série temporal de contagem
Velicer & Fava, 2003
“Em seguida, estudamos o número mensal de motoristas mortos de veículos leves na Grã-Bretanha entre janeiro de 1969 e dezembro de 1984 mostrado na Figura 3.
Esta série temporal é parte de um conjunto de dados que foi considerado pela primeira vez por Harvey e Durbin (1986) para estudar o efeito de uso obrigatório do cinto de segurança introduzido em 31 de janeiro de 1983.
O conjunto de dados, incluindo covariáveis adicionais, está disponível em R no objeto Seatbelts. Em seu artigo Harvey e Durbin (1986) analisam os números de acidentes para motoristas e passageiros de carros, que são tão grandes que podem ser tratados com métodos para dados de valor contínuo.
O número mensal de motoristas mortos de vans analisados aqui é muito menor (seu mínimo é 2 e seu máximo 17) e, portanto, métodos para dados de contagem devem ser preferidos.
Para a seleção do modelo, usamos apenas os dados até dezembro de 1981.”
Liboschik et al., 2017, p. 18
Liboschik et al., 2017
Time-Series [1:192] from 1969 to 1985: 12 6 12 8 10 13 11 6 10 16 ... Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1969 12 6 12 8 10 13 11 6 10 16 13 14
1970 14 6 8 11 7 13 13 11 11 14 16 14
1971 17 16 15 13 13 15 12 6 9 13 14 15
1972 14 3 12 13 12 8 8 15 8 5 17 14
1973 13 5 8 5 12 11 13 15 11 11 10 13
1974 8 6 8 14 12 14 13 9 4 13 6 15
1975 12 16 7 12 10 9 9 6 7 13 14 13
1976 14 11 11 10 4 8 9 10 10 5 13 12
1977 10 9 7 5 10 5 6 8 6 12 15 7
1978 14 4 10 8 7 11 3 5 11 10 10 7
1979 10 11 9 7 8 13 8 5 8 7 12 10
1980 7 4 10 4 8 8 7 10 8 14 8 9
1981 8 6 7 6 5 4 5 10 7 10 12 7
1982 4 5 6 4 4 8 8 3 7 12 2 7
1983 8 3 2 6 3 7 6 8 8 4 3 5
1984 5 3 4 3 6 6 7 5 7 7 4 7plot(timeseries,
xlab="Mes", ylab="Numero de obitos",
main="Monthly number of killed van drivers in Great Britain",
col="black")
Created Using changepoint version 2.3
Changepoint type : Change in mean
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 15.77249
Minimum Segment Length : 1
Maximum no. of cpts : Inf
Number of changepoints: 39 
Created Using changepoint version 2.3
Changepoint type : Change in variance
Method of analysis : PELT
Test Statistic : Normal
Type of penalty : MBIC with value, 15.77249
Minimum Segment Length : 2
Maximum no. of cpts : Inf
Changepoint Locations : 
Optimal 3-segment partition:
Call:
breakpoints.formula(formula = timeseries ~ 1)
Breakpoints at observation number:
88 156
Corresponding to breakdates:
1976(4) 1981(12) 
regressors <- cbind(PetrolPrice = Seatbelts[, c("PetrolPrice")],
linearTrend = seq(along = timeseries)/12)
timeseries_until1981 <- window(timeseries, end = 1981 + 11/12)
regressors_until1981 <- window(regressors, end = 1981 + 11/12)
seatbeltsfit <- tscount::tsglm(timeseries_until1981,
link = "log",
distr = "poisson",
model = list(past_obs = c(1, 12)),
xreg = regressors_until1981)
summary(seatbeltsfit)
Call:
tscount::tsglm(ts = timeseries_until1981, model = list(past_obs = c(1,
12)), xreg = regressors_until1981, link = "log", distr = "poisson")
Coefficients:
Estimate Std.Error CI(lower) CI(upper)
(Intercept) 1.8347 0.36793 1.1135 2.5558
beta_1 0.0866 0.08092 -0.0720 0.2452
beta_12 0.1535 0.08453 -0.0122 0.3191
PetrolPrice 0.7787 2.30368 -3.7364 5.2939
linearTrend -0.0303 0.00784 -0.0456 -0.0149
Standard errors and confidence intervals (level = 95 %) obtained
by normal approximation.
Link function: log
Distribution family: poisson
Number of coefficients: 5
Log-likelihood: -396.1849
AIC: 802.3698
BIC: 817.6191
QIC: 802.3698
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.8346719 0.3679315 4.9865 6.15e-07 ***
beta_1 0.0865713 0.0809250 1.0698 0.2847216
beta_12 0.1534817 0.0845250 1.8158 0.0693991 .
PetrolPrice 0.7787267 2.3036794 0.3380 0.7353360
linearTrend -0.0302632 0.0078409 -3.8597 0.0001135 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1timeseries_1982 <- window(timeseries, start = 1982, end = 1982 + 11/12)
regressors_1982 <- window(regressors, start = 1982, end = 1982 + 11/12)
predict(seatbeltsfit,
n.ahead = 12,
level = 0.9,
global = TRUE,
B = 2000,
newxreg = regressors_1982)$pred Jan Feb Mar Apr May Jun Jul Aug
1982 7.719763 7.437318 7.555106 7.406202 7.203925 6.995341 7.155932 7.862240
Sep Oct Nov Dec
1982 7.533187 7.858005 8.064586 7.478864seatbeltsfit_alldata <- tscount::tsglm(timeseries,
link = "log",
distr = "poisson",
model = list(past_obs = c(1, 12)),
xreg = regressors)
tscount::interv_test(seatbeltsfit_alldata,
tau = 170,
delta = 1,
est_interv = TRUE)
Score test on intervention(s) of given type at given time
Chisq-Statistic: 1.108797 on 1 degree(s) of freedom, p-value: 0.2923436
Fitted model with the specified intervention:
Call:
tsglm(ts = fit$ts, model = model_extended, xreg = xreg_extended,
link = fit$link, distr = fit$distr)
Coefficients:
(Intercept) beta_1 beta_12 PetrolPrice linearTrend interv_1
1.93298 0.08178 0.13943 0.41863 -0.03466 -0.21683 “With a p value of 0.29 the null hypothesis of no intervention cannot be rejected at a 5% significance level.”
Liboschik et al., 2017, p. 21
Sequência em R
[1] 1 2 3[1] 1 2 3[1] 1 2 3[1] -3 -2 -1 0 1 2 3[1] -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0[1] -10 -5 0 5 10[1] -10.0 -7.5 -5.0 -2.5 0.0 2.5 5.0 7.5 10.0[1] 0.0 0.2 0.4 0.6 0.8 1.0[1] 1 2 4 8 16 32# 1 2 4 8 16 32
geomSeries <- function(base, max){
base^(0:floor(log(max, base)))
}
geomSeries(base=2, max=50)[1] 1 2 4 8 16 32# 1 2 4 8 16 32
geomSeq <- function(start,ratio,begin,end){
begin<-begin-1
end<-end-1
start*ratio**(begin:end)
}
geomSeq(1, 2, 1, 6)[1] 1 2 4 8 16 32[1] 4 8 16 32[1] "2018-01-01" "2019-01-01" "2020-01-01" "2021-01-01"# "2018-01-01" "2019-01-01" "2020-01-01" "2021-01-01"
# month series
seq(as.Date("2000/1/1"), by="month", length=12) [1] "2000-01-01" "2000-02-01" "2000-03-01" "2000-04-01" "2000-05-01"
[6] "2000-06-01" "2000-07-01" "2000-08-01" "2000-09-01" "2000-10-01"
[11] "2000-11-01" "2000-12-01"# [1] "2000-01-01" "2000-02-01" "2000-03-01" "2000-04-01"
# [5] "2000-05-01" "2000-06-01" "2000-07-01" "2000-08-01"
# [9] "2000-09-01" "2000-10-01" "2000-11-01" "2000-12-01"
# quarter series
seq(as.Date("2000/1/1"), as.Date("2001/1/1"), by="3 months")[1] "2000-01-01" "2000-04-01" "2000-07-01" "2000-10-01" "2001-01-01"# [1] "2000-01-01" "2000-04-01" "2000-07-01" "2000-10-01"
# [5] "2001-01-01"
# week series
seq(as.Date("2001/1/1"), as.Date("2001/3/1"), by="weeks")[1] "2001-01-01" "2001-01-08" "2001-01-15" "2001-01-22" "2001-01-29"
[6] "2001-02-05" "2001-02-12" "2001-02-19" "2001-02-26"# [1] "2001-01-01" "2001-01-08" "2001-01-15" "2001-01-22"
# [5] "2001-01-29" "2001-02-05" "2001-02-12" "2001-02-19"
# [9] "2001-02-26"
# day series
seq(as.Date("2001/1/1"), as.Date("2001/1/31"), by="days") [1] "2001-01-01" "2001-01-02" "2001-01-03" "2001-01-04" "2001-01-05"
[6] "2001-01-06" "2001-01-07" "2001-01-08" "2001-01-09" "2001-01-10"
[11] "2001-01-11" "2001-01-12" "2001-01-13" "2001-01-14" "2001-01-15"
[16] "2001-01-16" "2001-01-17" "2001-01-18" "2001-01-19" "2001-01-20"
[21] "2001-01-21" "2001-01-22" "2001-01-23" "2001-01-24" "2001-01-25"
[26] "2001-01-26" "2001-01-27" "2001-01-28" "2001-01-29" "2001-01-30"
[31] "2001-01-31"# [1] "2001-01-01" "2001-01-02" "2001-01-03" "2001-01-04"
# [5] "2001-01-05" "2001-01-06" "2001-01-07" "2001-01-08"
# [9] "2001-01-09" "2001-01-10" "2001-01-11" "2001-01-12"
# [13] "2001-01-13" "2001-01-14" "2001-01-15" "2001-01-16"
# [17] "2001-01-17" "2001-01-18" "2001-01-19" "2001-01-20"
# [21] "2001-01-21" "2001-01-22" "2001-01-23" "2001-01-24"
# [25] "2001-01-25" "2001-01-26" "2001-01-27" "2001-01-28"
# [29] "2001-01-29" "2001-01-30" "2001-01-31"
# first four months
cat(month.name[1:4], sep=" ")January February March AprilJanuary-February-March-AprilJanuaryFebruaryMarchAprilStatistics in R [,1] [,2] [,3] [,4] [,5]
[1,] "Jan 1999" "Jan 2000" "Jan 2001" "Jan 2002" "Jan 2003"
[2,] "Feb 1999" "Feb 2000" "Feb 2001" "Feb 2002" "Feb 2003"
[3,] "Mar 1999" "Mar 2000" "Mar 2001" "Mar 2002" "Mar 2003"
[4,] "Apr 1999" "Apr 2000" "Apr 2001" "Apr 2002" "Apr 2003"
[5,] "May 1999" "May 2000" "May 2001" "May 2002" "May 2003"
[6,] "Jun 1999" "Jun 2000" "Jun 2001" "Jun 2002" "Jun 2003"
[7,] "Jul 1999" "Jul 2000" "Jul 2001" "Jul 2002" "Jul 2003"
[8,] "Aug 1999" "Aug 2000" "Aug 2001" "Aug 2002" "Aug 2003"
[9,] "Sep 1999" "Sep 2000" "Sep 2001" "Sep 2002" "Sep 2003"
[10,] "Oct 1999" "Oct 2000" "Oct 2001" "Oct 2002" "Oct 2003"
[11,] "Nov 1999" "Nov 2000" "Nov 2001" "Nov 2002" "Nov 2003"
[12,] "Dec 1999" "Dec 2000" "Dec 2001" "Dec 2002" "Dec 2003"[1] 1 1[1] 4 3 4 3[1] 4 4 3 3 3 3 3[1] 4 4 3 3[1] "a" "a"[1] "a" "b" "a" "b"[1] "a" "a" "b" "b" "b" "b" "b"[1] "a" "a" "b" "b"[1] "1" "2" "3" "4" "5" "6"[1] "1 st" "2 nd" "3 rd" "4 th" "5 th"[1] "1st" "2nd" "3rd" "4th" "5th"[1] "1st, 2nd, 3rd, 4th, 5th"[1] "Q1" "Q2" "Q3" "Q4" "Q5" "Q6"[1] "Q1" "Q2" "Q3" "Q4" "Q5" "Q6"[1] "Today is Fri Aug 8 20:41:12 2025"[1] "Today is Fri Aug 8 20:41:12 2025"[1] "Q1+Q2+Q3+Q4+Q5+Q6"[1] "F~Q1+Q2+Q3+Q4+Q5+Q6"[1] "F~Q1+Q2+Q3+Q4+Q5+Q6\nG~Q7+Q8+Q9"
(-6,-5] (-5,-4] (-4,-3] (-3,-2] (-2,-1] (-1,0] (0,1] (1,2] (2,3] (3,4]
1 29 1281 21433 135303 342086 341153 135832 21597 1256
(4,5] (5,6]
29 0
(-6,-5] (-5,-4] (-4,-3] (-3,-2] (-2,-1] (-1,0] (0,1] (1,2]
0.000001 0.000029 0.001281 0.021433 0.135303 0.342086 0.341153 0.135832
(2,3] (3,4] (4,5] (5,6]
0.021597 0.001256 0.000029 0.000000 [1] (0.994,3] (0.994,3] (0.994,3] (3,5] (3,5] (0.994,3] (0.994,3]
[8] (3,5] (3,5] (5,7.01] (5,7.01]
Levels: (0.994,3] (3,5] (5,7.01] [1] (0.994,3] (0.994,3] (0.994,3] (3,5] (3,5] (0.994,3] (0.994,3]
[8] (3,5] (3,5] (5,7.01] (5,7.01]
Levels: (0.994,3] < (3,5] < (5,7.01]## one way to extract the breakpoints
labs <- levels(cut(aaa, 3))
cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs) )) lower upper
[1,] 0.994 3.00
[2,] 3.000 5.00
[3,] 5.000 7.01[1] Control Control Control Control Treat Treat Treat Treat
Levels: Control Treat# [1] Control Control Control Control Treat Treat
# [7] Treat Treat
# Levels: Control Treat
## 20 alternating 1s and 2s
gl(n=2, k=1, length=2*4)[1] 1 2 1 2 1 2 1 2
Levels: 1 2[1] 1 1 2 2 1 1 2 2
Levels: 1 2 [1] Low Low Low Low Medium Medium Medium Medium High High
[11] High High
Levels: Low < Medium < HighReferências
Andrews, DF & Herzberg, AM (1985) Data: A collection of problems from many fields for the student and research worker. NY: Springer-Verlag.
ANTUNES, JLF e CARDOSO, MRAs. Uso da análise de séries temporais em estudos epidemiológicos. Epidemiol. Serv. Saúde [online]. 2015, vol.24, n.3 [citado 2023-02-12], pp.565-576. Disponível em: http://scielo.iec.gov.br/scielo.php?script=sci_arttext&pid=S1679-49742015000300024&lng=pt&nrm=iso. ISSN 1679-4974.
Box, GEP, Jenkins, GM, Reinsel, GC, and Ljung, GM (2015) Time series analysis: Forecasting and control. 5th ed. NJ: Wiley.
Campbell, W et al. (1997) The econometrics of financial markets. 2nd ed. NJ: Princeton.
Campbell, DT & Stanley, JC (1969) Delineamentos experimentais e quase-experimentais de pesquisa. São Paulo: EPU/EDUSP. Capítulo 7: O experimento de série temporal.
Cunha, GZ (2016) Modelos de séries temporais para dados de contagem. Dissertação de mestrado do PBE-UEM. http://pbe.uem.br/wp-content/uploads/2020/03/MODELOS-DE-S%C3%89RIES-TEMPORAIS-PARA-DADOS-DE-CONTAGEM-%E2%80%93-Guilherme-Zubatch-da-Cunha.pdf.
Enders, W (2015) Applied econometric time series. 4th ed. NJ: Wiley. http://time-series.net/.
Ferreira, PGC et al. (2015) Box & Jenkins com função de transferência usando o R: um estudo de série de Energia Natural Afluente (ENA) do subsistema sul. SBPO.
Furrukh, M (2013) Tobacco smoking and lung cancer: Perception-changing facts. Sultan Qaboos Univ Med J. 13(3):345-58. doi: 10.12816/0003255.
Hamaker, EL (2013) Why researchers should think “within-person”: A paradigmatic rationale. In Mehl MR &, Conner TS (Eds.) (2013) Handbook of research methods for studying daily life. NY: Guilford. p. 43-61.
Hamaker EL, Asparouhov T, Brose A, Schmiedek F, Muthén B (2018) At the frontiers of modeling intensive longitudinal data: Dynamic structural equation models for the affective measurements from the COGITO Study. Multivariate Behav Res 53(6):820-41. doi: 10.1080/00273171.2018.1446819.
Hamilton, JD (1994) Time series analysis. NJ: Princeton.
Hecht, SS (1999) Tobacco Smoke Carcinogens and Lung Cancer. JNCI: Journal of the National Cancer Institute 91(14): 1194–210. https://doi.org/10.1093/jnci/91.14.1194.
Houghton, A, Munster, EW and Viola, MV (1978). Increased incidence of malignant melanoma after peaks of sunspot activity. The Lancet 8: 759-60.
IBM SPSS Forecasting 22 (2014). USA: IBM SPSS.
Liboschik, T; Fokianos, K & Fried, R (2017). tscount: An R package for analysis of count time series following generalized linear models. Journal of Statistical Software, 82(5): 1–51. https://doi.org/10.18637/jss.v082.i05.
Mair, P (2018) Modern psychometrics with R. USA: Springer. Chapter 13: Modeling trajectories and time series; Chapter 14: Analysis of fMRI data.
Mattson ME, Pollack ES, Cullen JW. (1987) What are the odds that smoking will kill you? Am J Public Health. 77(4):425-31. doi: 10.2105/ajph.77.4.425. Erratum in: Am J Public Health 1987 Jul;77(7):818.
Morgenstern, H (1995) Ecologic studies in epidemiology: concepts, principles, and methods. Annu Rev Public Health 16:61-81. doi: 10.1146/annurev.pu.16.050195.000425.
Murteira, BJF et al. (1993) Análise de sucessões cronológicas. Lisboa: McGraw-Hill. Capítulo 9: Modelos de função de transferência e de intervenção.
Newbold, P & Bos, T (1990) Introductory business forecasting. South-Western.
Stadnitski, T (2020) Time series analyses with psychometric data. PLoS ONE 15(4): e0231785. https://doi.org/10.1371/journal.
Morley, S (2018) Single-case methods in clinical psychology: A practical guide. NY: Routledge. Chapter 5: Visual analysis if single-case data; Chapter 6: Statistical analysis of single-case data.
Schaffer AL, Dobbins TA, Pearson SA (2021) Interrupted time series analysis using autoregressive integrated moving average (ARIMA) models: a guide for evaluating large-scale health interventions. BMC Med Res Methodol 21(1):58. doi: 10.1186/s12874-021-01235-8.
Spata, AV (2005) Métodos de pesquisa: ciências do comportamento e diversidade humana. RJ: LTC. Capítulo 13: Delineamentos de sujeito único.
Schlesselman, JJ (2006) The emerging case-control study: Lung cancer in relation to tobacco smoking. Preventive Medicine 43(4): 251-255. https://doi.org/10.1016/j.ypmed.2006.07.014.
Tsay, RS (2010) Analysis of financial time series. 3rd ed. NJ: Wiley. * Tsay, RS (2012) An introduction to analysis of financial data with R. NJ: Wiley. https://faculty.chicagobooth.edu/ruey-s-tsay/research/an-introduction-to-analysis-of-financial-data-with-r.
Tsay, RS (2014) Multivariate time series analysis with R and financial applications. NJ: Wiley. https://faculty.chicagobooth.edu/ruey-s-tsay/research/multivariate-time-series-analysis-with-r-and-financial-applications.
Tsay, RS & Chen, R (2018) Nonlinear time series analysis. NJ: Wiley. https://faculty.chicagobooth.edu/ruey-s-tsay/research/nonlinear-time-series-analysis.
Valachovic, E & Zurbenko, I (2014) Skin cancer, irradiation, and sunspots: The solar cycle effect. BioMed Research International 2014. https://doi.org/10.1155/2014/538574.
Velicer, WF & Fava, JL (2003) Time series analysis. In Schinka, J & Velicer, WF (Eds.), Research Methods in Psychology (581-606). Volume 2, Handbook of Psychology (Weiner, IB Editor-inChief). NY: Wiley.
Valadkhani, A (2025) Inflation-driven instability in US sectoral betas. J Asset Manag. https://doi.org/10.1057/s41260-025-00413-3
Victor-Edema, UA & Essi, DI. A transfer function modelling of Nigerian current account (Net) and exchange rate. Int J Stat Appl Math 5(4): 177-185. DOI: 10.22271/maths.2020.v5.i4c.564.
Zhang YQ, Xu ZT, Yu QL (1993) The carcinogenic effect of solar activity with different intensities on embryos. J Radiat Res. 34(3):240-6. doi: 10.1269/jrr.34.240. Erratum in: J Radiat Res (Tokyo) 1993 Dec;34(4):419.
Hyndman, RJ. The ARIMAX model muddle. https://robjhyndman.com/hyndsight/arimax/
Ren, R (2020) All the confusion about ARIMA, ARIMAX, Transfer Function, Dynamic Regression Models. https://ruqinren.wordpress.com/2020/02/21/all-the-confusion-about-arima-arimax-transfer-function-dynamic-regression-models/
ARIMA Intervention Transfer Function: How to Visualize the Effect : https://stats.stackexchange.com/questions/108374/arima-intervention-transfer-function-how-to-visualize-the-effect.