Series de Tiempo Miguel de Carvalho DEPARTAMENTO DEPART AMENTO DE ESTADÍSTIC ESTADÍSTIC A PONTIFICIA PONTIFICIA UNIVERSIDAD CATÓLICA CATÓLICA DE CHILE
[email protected]
§1. INTRODUCCIÓN AL
MODELAMIENTO ESTADÍSTICO ESTADÍSTICO DE SERIES DE TIEMPO
Parte I Introducción y Motivación
Parte I Introducción y Motivación
Introducción Las series de tiempo son un tema de amplio interés en estadística con diversas aplicaciones práticas. Mientras que en muchos contextos se parte del supuesto de que los datos son una muestra aleatoria (y por consiguiente independientes) Y t t
i.i.d.
∼ F (y ), Y Y
t = 1, . . . , n ,
en nuestro contexto los datos estarán ordenados en el tiempo y por tanto son a menudo dependientes. dependencia temporal Es así importante identificar la estructura de dependencia que se oculta en los datos, para por ejemplo
utilizar eventos pasados para predecir eventos futuros, o extraer el comportamiento cíclico de una serie de tiempo.
Economia Estadounidense—Mercado Laboral US Department of Labor—Employment & Training Administration http://www.ows.doleta.gov
0 0 7 0 0 6
) s e r a 0 l l i 0 5 m ( s m 0 i a l 0 4 C l a i t i n 0 I 0 3 0 0 2
1970
1980
1990
2000
Tiempo (años) †
Nota: Datos desestacionalizados.
2010
Economia Chilena—Estadísticas Macro Banco Central de Chile
http://si3.bcentral.cl
0 2 1 8 0 1
1 1
5 1 1
0 1
) % ( o e l p m e s e D e d a s a T
6 0 1
9
C E C A M I
0 1 1 C P I
4 0 1
5 0 1
8
2 0 1
7
0 0 1 0 0 1
6
2 01 01 0
2 01 01 1
2 01 01 2
Tiempo (años)
†
2 01 01 3
2008
2010
2012
Tiempo (años)
20 00 09
20 01 10
20 01 11
20 01 12
20 01 13
Tiempo (años)
IMACEC: Indicador Mensual de Actividad Económica; IPC: Índice de Precios al Consumidor. Consumidor.
Índice de Precio Selectivo de Acciones (IPSA) Bolsa de Comercio de Santiago
http://www.bolsantiago.cl
0 0 0 5
2 . 0
1 . 0
0 0 0 4
A S P I
A S P I s o n r o t e R g o l −
0 0 0 3
0 0 0 2
0 . 0
1 . 0 −
2 . 0 −
3 . 0 −
0 0 0 1
2004
2006
2008 Tiempo (años)
2010
2012
2004
2006
2008 Tiempo (años)
2010
2012
São Paulo, Brasil—Datos de Polución Companhia de Tecnologia de Saneamento Ambiental http://www.cetesb.sp.gov.br
Kazajstán—Granjas de Producción de Grano Granja 1 ) . n o t 1 . 0 ( n ó i c c u d o r P
) . n o t
1 . 0 ( n ó i c c u d o r P
5 1 0 1 5
1 . 0 ( n ó i c c u d o r
1 98 0
2 00 0
) . n o t 1 . 0 ( n ó i c c u d o r P
0 2 5 1 0 1 5
1 96 0
1 98 0
2 00 0
0 2
Granja 4 ) . n o t 1 . 0 ( n ó i c c u d o r P
5 1 0 1
5
1 96 0
1 98 0
2 00 0
0 2 5 1 0 1
5
1 96 0
1 98 0
2 00 0
Tiempo
Tiempo
Tiempo
Granja 5
Granja 6
Granja 7
Granja 8
) . n o t
1 . 0 ( n ó i c c u d o r P
5 1 0 1 5
1 98 0
2 00 0
) . n o t
5 1
1 . 0 ( n ó i c c u d o r P
0 1
5
1 96 0
1 98 0
2 00 0
Tiempo
Tiempo
Granja 9
Granja 10 ) . n o t
5 1
1 . 0 ( n ó i c c u d o r
0 1
5
P
1 . 0 ( n ó i c c u d o r
0 1
5
1 98 0 Tiempo
2 00 0
) . n o t
0 2
1 . 0 ( n ó i c c u d o r P
5 1 0 1 5
1 96 0
) . n o t
5 1
P
1 96 0
Granja 3
Tiempo
0 2
1 96 0
) . n o t
) . n o t 1 . 0 ( n ó i c c u d o r P
0 2
1 96 0
Granja 2
1 98 0
2 00 0
1 98 0 Tiempo
2 00 0
5 1 0 1
5
1 96 0
1 98 0
2 00 0
Tiempo
Tiempo
Granja 11
Granja 12
0 2
) . n o t
5 1
1 . 0 ( n ó i c c u d o r
0 1 5
P
1 96 0
0 2
P
1 96 0
1 98 0 Tiempo
2 00 0
0 2 5 1 0 1 5 0
1 96 0
1 98 0 Tiempo
2 00 0
Bosque de Beatenberg, Suiza—Datos Atmosféricos Langfristige Waldökosystem-Forschung http://www.wsl.ch
0 4
)
0 3
C °
(
a i r a i D a i d e M a r u t a r e p m e T
0 2
0 1
0
0 1 −
0 2 −
1998
2000
2002 Tiempo (años)
2004
2006
Alpes Suizos—Suelo Permafrost WSL Institute for Snow and Avalanche Research http://www.slf.ch
2
0 ) C º (
a r u 2 t a r − e p m e T 4 −
6 −
1998
2000
2002
2004
2006
2008
2010
Tiempo (años)
Figura : Temperatura del suelo sobre diferentes profundidades.
Salud Pública en Portugal—Influenza y Gripe Instituto Nacional de Saúde, Dr. Ricardo Jorge http://www.insa.pt
0 5 2 a n a 0 m 0 e 2 S r o p s 0 e 5 t r 1 e u M e 0 d 0 1 o r e m ú N 0 5
1980
1985
1990
1995
Tiempo (años)
2000
2005
Detalles sobre el Funcionamiento del Curso ∼
http://www.mat.puc.cl/ mdecarvalho/series_tiempo.html
Horario: martes y jueves, 17:00–18:20, sala 3. Evaluación:
examen final (60 %) proyecto de análisis de datos (40 %) mini-proyectos (bónus). Mini-proyectos. Los mini-proyectos serán cinco y funcionan como un
bonus —su entrega no es obligatoria. Proyecto de análisis de datos. Plazos:
Propuesta: 30 abril, 2013; Entrega de Proyecto: 4 junio, 2013. Referencias clave. Brockwell y Davis (2001), Shumway y Stoffer
(2011).
Parte II Introducción al Modelamiento Estadístico de Series de Tiempo —Dominio del Tiempo— El Análisis de la Dinamica Temporal de los Datos
Procesos Estocásticos Proceso Estocástico y Trayectoria
Definición (Proceso Estocástico) Sea T ⊆ R un conjunto de indexación. Un proceso estocástico es una familia de variables aleatorias {Y t }t T definidas en un espacio de probabilidad (Ω, A, P ). Para ω ∈ Ω fijo, la función y t ≡ Y t (ω), para t ∈ T, es una trayectoria del proceso estocástico {X t }t T . ∈
∈
Para series de tiempo igualmente espaciadas en el tiempo, el conjunto de indexación natural es T = Z, pero para otro tipo de datos T puede asumir otras formas (por ejemplo, T = R para datos continuos).
Características de un Proceso Estocástico Funciones de Media, Varianza, Autocovarianza y Autocorrelación
Definición Sea {Y t }t T un proceso estocástico dónde E |Y t | < ∞; además, sea F t (y ) = P (Y t ≤ y ). Las funciones de media y varianza son definida como ∈
µt = E (Y t ) =
y dF t (y ),
σt 2 = var(Y t ) = E (Y t
R
− µ )2 = t
−
(y µt )2 dF t (y ).
R
Definición (Funciones de Autocovarianza y Autocorrelación) Sea {Y t }t T un proceso estocástico donde σ t 2 < ∞, para t ∈ T. La función de autocovarianza γ (· , ·) de {X t }t T es definida como ∈
∈
γ (r , s ) = cov(Y r , Y s ) = E (Y r
{ − µ )(Y − µ )}, r , s ∈ T. La función de autocorrelación ρ (· , ·) de {X } es definida como Y
ρ(r , s ) = corr(Y r , Y s ) =
{
r
s
s
t t ∈T
γ (r , s ) γ (s , s )γ (t , t )
}1 2 , /
r , s
∈ T.
Estacionariedad
Estacionariedad—Definiciones
Definición Una serie de tiempo {Y t }t Z es fuertemente estacionaria si para todo el h ∈ Z y t 1 , . . . , t k ∈ Z ∈
d
(Y t1 , . . . , Y t k ) = (Y t1 +h , . . . , Y tk +h ).
Una serie de tiempo {Y t }t Z es (débilmente) estacionaria si, ∈
1 E Y t 2 <
∞, para todo el t ∈ Z, = µ, para todo el t ∈ Z,
| |
2 µt
3 γ (r , s ) = γ (r + t , s + t ), para todo el r , s , t
∈ Z.
Estacionariedad
Estacionariedad—Intuición y Otras Consideraciones En palavras: Una serie de tiempo es fuertemente estacionaria si las distribuciones conjuntas son invariantes a cambios de tiempo. Una serie de tiempo es (débilmente) estacionaria si fluctúa alrededor de un nivel fijo (µ) y la dependencia de observaciones pasadas depende sólo de su rezago. Si {Y t }t Z es estacionario podemos redefinir la autocovarianza y autocorrelación del proceso como funciones de una única variable h ∈ Z, ∈
≡ γ (0, h ) = cov(X , X
γ h
t
t +h ),
≡ ρ(0, h ) = corr(X , X
ρh
t
t +h ).
En la práctica, a menudo nos enfrentamos con la cuestión de si los datos son una trayectoria de un proceso estacionario; esto puede ser estudiado usando R y la prueba de KPSS, que será introducida adelante.
Ruido Blanco Definición
Definición Un proceso estocástico {Y t } es llamado ruido blanco si ninguno de sus elementos es correlacionado, y si µt = 0, para todo el t ∈ Z, σt 2 = σ 2 , para todo el t ∈ Z. Si además los Y t siguen una distribución normal, entonces se denomina proceso de ruido blanco Gaussiano, Y t
i.i.d.
∼ N (0, σ2 ).
El término ‘blanco’ viene de la analogía con la luz blanca, y indica que todas las frecuencias están igualmente presentes.
Ejemplo Simulado
Ruido Blanco Gaussiano Estándar ## set.seed es usado para efectos de reproducibilidad set.seed(100) n <- 1000 y <- rnorm(n) plot(y) Ruido Blanco, Yt ~ N(0,1)
3
2
1
t
0
1 − 2 − 3 −
0
200
400
600
Tiempo (t)
800
1000
Funciones de Autocovarianza Propriedades Teóricas
Teorema 1 Si γ ( ) es una función de autocovarianza de un proceso estocástico
{Y }
·
t t ∈Z ,
∈Z |γ | ≤ γ 0,
entonces para todo h
≥ 0, γ = γ . Una función g : Z → R es una función de autocovarianza de una serie de tiempo estacionaria si y sólo si es par (g (h ) = g (−h ) ) y no-negativa γ 0
2
definida , a saber
h
n
−h
n
− t )a ≥ 0,
a i g (t i
i =1 j =1
para todo n
h
T
j
n
j
T
n
∈ N, a = (a 1 , . . . , a ) ∈ R , y t = (t 1, . . . , t ) ∈ Z . n
n
Función de Autocorrelación Estimación
Ya que ρ h = ρ h en la prática consideramos sólo las estimaciones de la función de autocorrelación de h = 0 hasta un determinado H ∈ N. −
El resultado es un gráfico que llamamos correlograma
{(h , ρ ) : h ∈ {0, . . . , H }},
h
dónde ρh es la función de autocorrelación muestral, construida usando los datos ( y 1 , . . . , y n ),
− − ρh =
n −h y )(y t +h t =1 ( y t n y )2 t =1 (y t
− y ) ,
1
−
y = n
n
t =1
y i .
Si y 1 , . . . , y n son realizaciones de un proceso de ruido blanco, con E (Y t 4 ) < ∞, entonces cuando n → ∞
√ n ρ ∼ N (0, 1).
·
h
Ver Brockwell y Davis (2001, Teorema 7.2.2). Para datos con tendencia, |ρh | va a exhibir un decaimiento lento; para datos con una componente periódica notable, ρh va a imitar el comportamiento cíclico con la misma periodicidad
Función de Autocorrelación R: acf{stats}
√
Las líneas azules corresponden a ±1,96/ n y son justificadas por el resultado asintótico de la diapositiva anterior. datos <- read.csv("ipsa.csv", header = TRUE) ## p.cierre: Precios de Cierre Mensuales, Bolsa de Mercado de Santiago p.cierre <- ts(rev(datos$Close), start = c(2003, 10), frequency = 52) acf(p.cierre, lag.max = 200) ## Bandas de confianza construidas de acuerdo con Brockwell and Davis (2001, Teorema 7.2.2). 0 . 1
8 . n 0 ó i c a l e 6 r . r 0 o c o t u A 4 . e 0 d n ó i c 2 n . u 0 F 0 . 0
0
1
2
3
4
Ejemplos
Paseo Aleatorio con Drift
Ejemplo 1 Sea T =
N0
. Sea ε t ∼ (0, σ 2 ) ruido blanco, Y 0 = 0 y defina Y t = µ + Y t −1 + εt ,
∈ N.
t
Muestre que {Y t } no es una serie de tiempo estacionaria. 2 Simule en R una trayectoria de {Y t }, con drift µ = 1. Paseo Aleatorio con Drift ( µ = 1)
0 0 8
0 0 6 t
0 0 4
0 0 2
0
0
200
400
600
Tiempo (t)
800
1000
Ejemplos
Proceso Alterno
Ejemplo Sea {εt } estacionario y Y t =
εt , t par, α + βεt , t impar.
con (α, β ) ∈ R2 . Verifique si existe (α, β ) ∈ R2 tal que {Y t } sea estacionario.
Ejemplos
Proceso Sinusoidal
Ejemplo 1 Sea Y t = Π1 cos(θt ) + Π2 sin(θt ),
θ
con Π 1 (i) Π2 y d
∼ (µ , σ2 ), , σ2 ) ∈ R × R
Π1 = Π2
Π
Π
(µΠ , σΠ2 )
∈ [−π, π], +
∈R×R . , tal que {Y } sea estacionario.
+ Verifique si existe (µΠ Π t 2 Use el ejercicio anterior para obtener un ejemplo de una función no-negativa definida.
Ejemplos
Procesos Lineales
Ejemplo 1 Defina un proceso lineal Y t como un proceso que puede ser escrito como una combinación lineal de ruido blanco ε t (0, σ 2 ),
∼
∞
Y t = µ +
| | ∞
Ψ j εt − j ,
j =−∞
Ψ j <
j =−∞
∞.
Muestre que
∞
2
γ h = σε
Ψ j Ψ j +h .
j =−∞
2 Muestre que los siguientes procesos son lineales y obtenga su
representación lineal Ruido Blanco: Y t = ε t , εt ∼ (0, σ 2 ). AR(1): Y t = β Y t 1 + εt , εt ∼ (0, σ 2 ), con |β | = ±1. = 0. MA(1): Y t = θεt 1 + εt , con ε t ∼ (0, σ2 ) y θ −
−
Ejemplos Simulados
Ruido Blanco Gaussiano Estándar
Ruido Blanco, Y t ~ N(0,1) 0 . 1
3
8 . n 0 ó i c a l e 6 r . r o 0 c o t u A 4 e . d 0 n ó i c n 2 u . F 0
2
1
t
0
1 − 2 −
0 . 0
3 −
0
200
400
600
Tiempo (t)
800
1000
0
20
40
60
Rezago
80
100
Ejemplos Simulados
Proceso Autoregresivo: AR(1), β = 0,1
4
0 . 1
3 8 . n 0 ó i c a l e 6 r . r 0 o c o t u A 4 . e 0 d n ó i c n 2 . u 0 F
2
1 t
0
1 − 2 −
0 . 0
3 −
0
200
400
600
Tiempo (t)
800
1000
0
20
40
60
Rezago
80
100
Ejemplos Simulados
Proceso Autoregresivo: AR(1), β = 0,9
0 . 1
8 . 0
n ó i c a 6 l . e 0 r r o c o t 4 u . A 0 e d n 2 . ó i 0 c n u F 0 . 0
5
t
0
5 −
2 . 0 −
0
200
400
600
Tiempo (t)
800
1000
0
20
40
60
Rezago
80
100
Ejemplos Simulados
Proceso Autoregresivo: MA(1), θ = 0,1
4
0 . 1
3 8 . n 0 ó i c a l e 6 r . r 0 o c o t u A 4 . e 0 d n ó i c n 2 . u 0 F
2
1 t
0
1 − 2 −
0 . 0
3 −
0
200
400
600
Tiempo (t)
800
1000
0
20
40
60
Rezago
80
100
Ejemplos Simulados
Proceso Autoregresivo: MA(1), θ = 0,9
0 . 1
4 8 . n 0 ó i c a l e 6 r . r o 0 c o t u A 4 . e 0 d n ó i c n 2 . u 0 F
2
t
0
2 −
0 . 0
4 −
0
200
400
600
Tiempo (t)
800
1000
0
20
40
60
Rezago
80
100
Inferencia sobre la Significancia Global de los Rezagos Para probar la significancia global de los primeros m rezagos H0 : ρ 1 , . . . , ρH = 0
= 0, Ha : ∃ i ∈ {1, . . . , H } : ρ i
vs .
podemos utilizar la prueba de Portmanteau (Box & Pierce, 1970) P = n Q H
H
∼ − ρ2h
·
χ2H ,
h =1
o la prueba de Ljung–Box (Ljung & Box, 1978) LB
H
Q H = n (n + 2)
h =1
ρh 2 n h
∼ χ2 . ·
H
El valor de H es seleccionado de un modo arbitrario; usualmente H = 20.
Prueba de Ruido Blanco Box.test{stats} + EuStockMarkets{datasets}
En finanzas es de interés probar la hipótesis de ruido blanco ya que corresponde a la célebre hipótesis de mercado eficiente (Fama, 1970; Lo and MacKinlay, 1999). Vamos ilustrar las ideas con algunos datos del Swiss Market Index (SMI), el principal índice de la bolsa de Suiza; véase http://www.six-swiss-exchange.com
El código R es lo siguiente: ## Daily Closing Prices of Major European Stock Indices, 1991--1998 ## Germany DAX (Ibis), Switzerland SMI, France CAC, and UK FTSE data(EuStockMarkets) SMI <- diff(log((EuStockMarkets[,2]))) # Box.test(SMI, lag = 10, type = "Box-Pierce") > Box.test(SMI, lag = 10, type = "Ljung-Box") Box-Ljung test data: SMI X-squared = 12.4887, df = 10, p-value = 0.2537
Prueba de Ruido Blanco Box.test{stats} + EuStockMarkets{datasets}
valores.p <- as.vector(NA) m <- 20 for(i in 1:m) { valores.p[i] <- Box.test(SMI, lag = i, type = "Ljung-Box")$p.value } plot(valores.p, ylab = "p-values", xlab = "Rezago") point(valores.p[10], pch = 10) 0 . 1 x o B . . . 8 g . n 0 u j L e d 6 a . b 0 e u r p a 4 l . a r 0 a p p − s 2 . e r 0 o l a V 0 . 0
5
10
15
20
Prueba de Normalidad R: jarque.bera.test{tseries}
Para probar la hipótesis de normalidad podemos usar por ejemplo la prueba de Jarque–Bera cuyo estadístico es
∼ − − JB = γ
n
6
2
+
κ
2
n
·
24
χ22 ,
y donde γ y κ representan la asimetría y (exceso de) curtosis muestrales 1
n
γ =
1
n
n t =1 (Y t
n t =1 (Y t
Y )3
Y )2
1
, 3/2
κ =
n
1
n
− Y )4 − 3. 2 2 1 (Y − Y )
n t =1 (Y t
n t =
t
Prueba de Estacionariedad R: kpss.test{tseries}
La prueba de KPSS ha sido propuesto por Kwiatkowski et al. (1992). Se basa en un modelo simple que descompone la serie en una tendencia determinista, un paseo aleatorio, y una serie de ruido blanco Y t = β t + R t + εt ,
R t = R t −1 + U t ,
U t
∼iid (0, σ2 ), U
εt
∼iid (0, σ2). ε
El interés está en H0 : σU 2 = 0, y dependiendo de si fijamos β = 0 o no, la hipótesis nula es de estacionariedad en nivel o estacionariedad en tendencia, respectivamente. El estadístico de la prueba es n
n −2 S t 2 KPSSh = 2 , σ h t =1
t
S t =
e i ,
t = 1, . . . , n .
i =1
Los e 1 , . . . , e n son los residuos de la regresión y t = α + I (β = 0)t + εt , y σh 2 es un estimador consistente de la varianza de largo plazo σ 2 = limn n 1 E (S n 2 ), basado en los residuos truncados en el rezago
h
−
→∞
Prueba de Estacionariedad R: kpss.test{tseries}
En kpss.test, el estimador consistente de σ 2 es el estimador de √ Newey–West (Newey and West, 1987) y por defecto h n = (3 n /13), donde · es la función piso. La distribución asintótica de KPSSh bajo la hipótesis nula es complicada y por ejemplo en el caso de estacionariedad en nivel
n
n −2 S t 2 KPSSh = σh 2 t =1
1
V 2 (r )dr ,
h n = o (n −1/2 ),
0
donde V (r ) ≡ B (r ) − rB (1) es un puente Browniano estándar, donde B (·) denota el movimiento Browniano estándar. La distribución bajo la hipótesis nula de estacionariedad tendencia es aún más complicada —es una función del llamado puente Browniano de segundo nivel—, y por lo tanto en la práctica los valores críticos deben ser obtenidos por simulación de Monte Carlo.
Prueba de Estacionariedad R: kpss.test{tseries}
## install.packages("tis") ## install.packages("tseries") require(tis) require(tseries) claims <- read.csv2("claims.csv") Claims <- tis(claims[,3], start = ti(19670107, tif = "wsaturday")) ; nc <-length(c) kpss.test(Claims) > kpss.test(Claims) KPSS Test for Level Stationarity data: claims$S.A. KPSS Level = 2.0924, Truncation lag parameter = 10, p-value = 0.01 Warning message: In kpss.test(claims$S.A.) : p-value smaller than printed p-value
H0 es rechazada: los datos no son estacionarios en nivel.
Parte III Introducción al Modelamiento Estadístico de Series de Tiempo Dominio de la Frecuencia —El Análisis de Movimientos Regulares—
Análisis Espectral Introducción
Vamos ahora abandonar el dominio del tiempo , y hacer una breve incursión en el dominio de la frecuencia , que es particularmente importante para el análisis de ciclos y movimientos regulares en los datos. Los movimientos de alta frecuencia son asociados a oscilaciones más rápidas mientras que los movimientos de baja frecuencia son asociados a oscilaciones más lentas. Para describir este tipo de fenómenos necesitamos utilizar modelos con una componente sinusoidal y recurrir a conceptos de análisis de Fourier. Vamos primero empezar por comprender el campo de aplicación de los métodos.
¿Porqué Debo Aprender Análisis Espectral?
Para ayudar a los economistas y financieros a entender la dinámica de las contracciones económicas
P T
4 0 . 0
P T
P T
P
T
P
T
P T P
T
P T
P T
P
T
2 0 . 0
0 0 . 0
2 0 . 0 -
4 0 . 0 -
6 0 . 0 -
1950
1960
1970
1980
1990
2000
Time
Figura : Portada de la The Economist en Abril 2011 (izquierda); análisis espectral del ciclo económico de la economia estadounidense (derecha).
2010
¿Porqué Debo Aprender Análisis Espectral? Para ayudar a los Médicos a entender la recurrencia de las pandemias
0 5 2 a n a m e S r o p s e t r e u M e d o r e m ú N
0 0 2
0 5 1
0 0 1
0 5
1980
1985
1990
1995
2000
Tiempo (años)
Figura : Portada de la Time en 2009 (izquierda); número de muertes por semana en Portugal debidas a Influenza y Gripe (derecha).
2005
¿Porque Debo Aprender Análisis Espectral?
Para ayudar a los Geólogos a analizar los movimientos periódicos en el suelo permafrost
2
0 ) C º (
a r u 2 t a r − e p m e T 4 −
6 −
2002
2004
2006
2008
2010
Tiempo (años)
Figura : Temperatura en suelo permafrost medida 2 metros bajo del suelo.
Densidad Espectral Introducción
Comenzamos con la introducción de la densidad espectral de un proceso estacionario que se define como
∞
f (ω) =
h =−∞
γ h e−2h πω i ,
∈ [−1/2; 1/2].
ω
La definición en el intervalo [ −1/2; 1/2] es pura conveniencia, y se puede también definir f en otros compactos. Como f (ω) = f (−ω) en la prática apenas representamos f (ω) en ω ∈ [0; 1/2]. Al igual que en el dominio del tiempo —donde la autocorrelación fue la función ‘natural’ para describir la dinámica de la serie al largo del tiempo— la densidad espectral es la función ‘natural’ para describir las características de recurrencia de una serie de tiempo.
Punto de Partida Teórico Teorema de Herglotz
Teorema (Herglotz) Una función γ h , h como
∈ Z, es no-negativa definida si y sólo si puede ser escrita
1/2
γ h =
2h πω i
e
dF (ω),
1/2
−
donde F : [ 1/2, 1/2] R es continua a la derecha y limitada, y unicamente determinada por las condiciones
−
−
1 F ( 1/2) = 0;
−
2 F (1/2) = γ 0 .
→
Punto de Partida Teórico
Elo entre la Densidad Espectral y la Función Autocovarianza de una Serie Estacionaria
Teorema Si γ h es la función de autocovarianza de un proceso estacionario Y t con h =∞
| |
γ h <
h =−∞
∞,
entonces la densidad espectral de Y t , es dada por dF dω
∞
(ω)
≡ f (ω) =
γ h e−2h πω i
h =−∞
= γ 0 +
∞
h =1
γ h cos(2h πω),
∈ [−1/2; 1/2].
ω
Periodograma
Estimación de la Densidad Espectral Defina la ordenada del periodograma en ω ∈ [0, 1/2) para una serie de tiempo igualmente espaciada I(ω)
=
1
n
2
n
Y t sin(2t πω)
t =1
+
2
n
Y t cos(2t πω)
t =1
El periodograma es un estimador de la densidad espectral muy frecuentemente utilizado, que se define el logaritmo del periodograma ordenadas evaluado en las frecuencias de Fourier, es decir
{(2 j π/n , log I( j /n )) : j = 1 + I (n impar), . . . , m = (n − 1)/2}. El periodograma acumulativo permite probar la hipótesis de ruido blanco, a través de la comparación de C r =
r j =1 I( j /n ) , m I ( / ) j n l =1
r = 1, . . . , m ,
con una distribución uniforme. La intuición es la siguiente: El ruido blanco tendrá una densidad espectral plana de manera que su periodograma acumulativo corresponde a la línea diagonal.
Computing the Spectrum in R Con Nuestro Código
## clear variables in memory rm(list = ls(all = T)) set.seed(1) y <- arima.sim(n = 99, list(ar = c(0.9,0))) y <- ts(y, frequency = 2 *pi) n <- length(y) ## ordenada del periodograma I <- function(w) { time <- seq(1, n) return(1/n*((sum(y*sin(2*pi*w*time)))^2 + (sum(y*cos(2*pi*w*time)))^2)) } ## Fourier frequencies m <- ceiling((n-1)/2) J <- seq(1, m) ff <- 2*pi/n*J I.jn <- rep(NA, m) log.I <- rep(NA, m) for(j in 1:m){ I.jn[j] <- I(j/n) log.I[j] <- log(I(j/n)) } par(mfrow = c(1,2)) plot(ff, I.jn/10, type = "l", xlab = "Frecuencias de Fourier", ylab ession(paste(log," ",italic(I(j/n))),")"), log
"y")
Estimación del Espectro en R Con Nuestro Código
0 0 + e 5
1 0 − e 5 )
n j I
(
g o l
2 0 − e 5
3 0 − e 5
0.0
0.5
1.0
1.5
2.0
Frecuencias de Fourier
2.5
3.0
Estimación del Espectro en R R: spec.pgram{stats}
Vamos ahora a comparar el resultado con la función de R spec.pgram(y, fast = FALSE, taper = 0, detrend = FALSE) Series: y Raw Periodogram
0 0 + e 5
m u r t c e p s
1 0 − e 5
2 0 − e 5
3 0 − e 5
0.0
0.5
1.0
1.5
2.0
frequency bandwidth 0.0183
2.5
3.0
Estimación del Espectro en R R: cpgram{stats}
Para calcular el periodograma acumulativo con nuestro código basta escribir cr <- cumsum(I.jn)/sum(I.jn)
Vamos ahora comparar con la función cpgram{stats} ya implementada en R par(mfrow = c(1,2), pty = "s") plot(c(0,ff), c(0,cr), type ="S", main = "Nuestro Periodograma Acumulativo", xlab = "Frecuencias de Fourier", ylab = "", xlim = c(0,pi), ylim = c(0,1)) cpgram(y, taper = 0)
Con Nuestro Código 0 . 1
0 . 1 o v i t a l u m u c A o r t c e p s E
Usando cpgram de R
8 . 0
8 . 0
6 . 0
6 . 0
4 . 0
4 . 0
2 . 0
2 . 0
0 . 0
0 . 0
0.0
1.0
2.0
3.0
0.0
1.0
2.0
3.0
Ejemplo Simulado
Ruido Blanco Gaussiano Estándar
Ruido Blanco, Y t ~ N(0,1)
1 0 + e 1
3
0 0 + e 1
2
1
t
1 0 o − r t e c 1 e p s E 2 0 − e 1
0
1 − 2 −
3 0 − e 1
3 −
0
200
400
600
Tiempo (t)
800
1000
0.0
0.1
0.2
0.3
Frecuencia bandwidth = 0.000289
0.4
0.5
Ejemplo Simulado Ruido Blanco
0 . 1
8 . 0
o v i t a l u m 6 . u 0 c A a m a r 4 . g 0 o d o i r e P 2 . 0
0 . 0 0.0
0.1
0.2
0.3
Frecuencia
0.4
0.5
Ejemplos en Datos Reales
Economia Estadounidense—Ciclo Económico
7 0 + e 1
5 0 + e 1
0 o 0 t u 0 r 0 B 1 o n 0 r 0 e t 0 n 8 I o t c 0 u 0 d 0 o 6 r P 0 0 0 4
o r t c e p 3 s 0 E + e 1
1 0 + e 1
0 0 0 2
1950
1960
1970
1980
1990
Tiempo (años)
2000
2010
0.0
0.5
1.0 Frecuencia bandwidth = 0.00481
1.5
2.0
Ejemplos en Datos Reales Alpes Suizos—Suelo Permafrost
0 5 2
3 0 + e 1
a n a 0 m 0 e 2 S r o p s 0 e 5 t r 1 e u M e 0 d 0 1 o r e m ú N 0 5
1 0 + e o 1 r t c e p s E 1 0 − e 1
3 0 − e 1
1980
1985
1990
1995
Tiempo (años)
2000
2005
0
5
10
15
Frecuencia bandwidth = 0.0117
20
25
Ejemplos en Datos Reales Alpes Suizos—Suelo Permafrost
2 2 0 + e 1
0 ) C º (
0 0 + e o 1 r t c e 2 p 0 s − E e 1
a r u 2 t a r − e p m e T 4 −
4 0 − e 1 6 0 − e 1
6 −
2002
2004
2006
Tiempo (años)
2008
2010
0.0
0.1
0.2
0.3
Frecuencia bandwidth = 7.7e−05
0.4
0.5
Parte IV Resumen