Análisis espectral con R

Francisco Parra Rodríguez
Función gdf(a)
Transforma los datos del vector a del dominio del tiempo al dominio de la frecuencia pre-multiplicandolos por la matriz ortogonal, A, sugerida por Harvey (1978)

gdf <- function(a) {
a <- matrix(a, nrow = 1)
n <- length(a)
uno <- as.numeric(1:n)
A <- matrix(rep(sqrt(1/n), n), nrow = 1)
for (i in 3:n - 1) {
if (i%%2 == 0) {
A1 <- matrix(sqrt(2/n) * cos(pi * i * (uno - 1)/n), nrow = 1)
A <- rbind(A, A1)
} else {
A2 <- matrix(sqrt(2/n) * sin(pi * (i - 1) * (uno - 1)/n), nrow = 1)
A <- rbind(A, A2)
}
}
AN <- matrix(sqrt(1/n) * (-1)^(uno + 1), nrow = 1)
A <- rbind(A, AN)
A %*% t(a)
}

Función periodograma (a)
Calcula y presenta el espectro de la serie “a”
periodograma <- function(a) {
cf <- gdf(a)
n <- length(a)
if (n%%2 == 0) {
m1 <- c(0)
m2 <- c()
for (i in 1:n) {
if (i%%2 == 0)
m1 <- c(m1, cf[i]) else m2 <- c(m2, cf[i])
}
m2 <- c(m2, 0)
frecuencia <- seq(0:(n/2))
frecuencia <- frecuencia - 1
omega <- pi * frecuencia/(n/2)
periodos <- n/frecuencia
densidad <- (m1^2 + m2^2)/(4 * pi)
tabla <- data.frame(omega, frecuencia, periodos, densidad)
tabla$densidad[(n/2 + 1)] <- 2 * tabla$densidad[(n/2 + 1)]
data.frame(tabla[2:(n/2 + 1), ])
} else {
m1 <- c(0)
m2 <- c()
for (i in 1:(n - 1)) {
if (i%%2 == 0)
m1 <- c(m1, cf[i]) else m2 <- c(m2, cf[i])
}
m2 <- c(m2, cf[n])
frecuencia <- seq(0:((n - 1)/2))
frecuencia <- frecuencia - 1
omega <- pi * frecuencia/(n/2)
periodos <- n/frecuencia
densidad <- (m1^2 + m2^2)/(4 * pi)
tabla <- data.frame(omega, frecuencia, periodos, densidad)
data.frame(tabla[2:((n + 1)/2), ])
}
}

Función gperiodogrma (a)
Presenta gráficamente el espectro de la variable “a”
gperiodograma <- function(a) {
tabla <- periodograma(a)
plot(tabla$frecuencia, tabla$densidad, main = "Espectro", ylab = "densidad",
xlab = "frecuencia", type = "l", col = "#ff0000")
}

Función td (a,b)
Realiza una prueba estadística para estudiar la dependencia serial sobre el periodograma acumulado de a, con una significación de 0,1(b=1); 0,05(b=2); 0,025(b=3); 0,01(b=4) y 0,005 (b=5) (Durbin; 1969)

td <- function(a, b) {
per <- periodograma(a)
p <- as.numeric(per$densidad)
n <- length(p)
s <- p[1]
t <- 1:n
for (i in 2:n) {
s1 <- p[i] + s[(i - 1)]
s <- c(s, s1)
s2 <- s/s[n]
}
while (n > 75) n <- 75
if (b == 1)
c <- Test[n, 1] else {
if (b == 2)
c <- Test[n, 2] else {
if (b == 2)
c <- Test[n, 3] else c <- Test[n, 4]
}
}
min <- -c + (t/length(p))
max <- c + (t/length(p))
data.frame(s2, min, max)
}

Necesita leer el fichero TD guardado como csv a partir de test de durbin

Test <- read.csv("TD.csv", header = TRUE, sep = ";", dec = ",")


Fuction gtd (a,b)
Presenta graficamente los resultados de la prueba de Durbin (Durbin; 1969) :

gtd <- function(a, b) {
S <- td(a, b)
plot(ts(S), plot.type = "single", lty = 1:3)
}

Enlace a RPubs:

Análisis espectral con R

test de durbin

Anuncios

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s