Hola Jorge:
Disculpa la tardanza pero me han tenido liado en otros menesteres. Yo no
soy ni muchísimo menos un experto ni en series temporales ni en R de hecho
retomo el tema después de muchos (demasiados) años aparcado y dedicándome
a labores de programación pura y dura.
Respecto a la diferencia de resultados con R y Statgraphics, no conozco el
proceso de selección del modelo ARIMA que hace Statgraphics por lo que no
puedo compararlo con el de R. Visto la experiencia que ha comentado
Gregorio en el hilo de los Arima automatizados, no parece que el
auto.arima de R sea un algoritmo del todo afinado, aun.
La diferencia básica entre uno y otro resultado, tal y como yo lo
interpreto, es que R aplica una diferencia regular (señal de que detecta
tendencia) y Statgraphics no.
En el gráfico que envías parece apreciarse una ligera tendencia creciente
(que habría que confirmar observando la acf u otros métodos) de la serie
para establecer si conviene o no realizar esa diferencia. Eso es el primer
punto de comprobación, que junto con el chequeo de homocedasticidad y
posible presencia de estacionalidad son los primeros pasos para hacer la
serie estacional y posteriormente determinar el modelo Arima.
Quiero decir que, una vez que hemos hecho la serie estacionaria (sin
tendencia, con la estacionalidad eliminada y con homocedasticidad), ya
entraríamos a determinar que ordenes p,q,P,Q son los que determinan el
mejor modelo.Aquí estriba la diferencia, R considera que no es
estacionaria y por eso aplica una diferencia regular para eliminar la
tendencia, mientras que Statgraphics si la considera estacionaria y
calcula un ARMA (1,0).
Todo esto es algo muy básico en el tratamiento de series temporales, que
seguro que conoces mas que de sobra, pero me ha parecido adecuado ponerlo
expresamente para tener claro el punto de partida.
Respecto a trabajos o documentos que te puedan ayudar no entiendo muy bien
a que te refieres con lo de "que sirva de referencia para conocer como
explicar la estimación de una serie", si me lo aclaras intentare
comentarte los documentos que yo manejo.
Por cierto el código de R que uso para conectarme al servidor de SQL,
bajarme los datos, calcular los arimas y guardar los resultados es mas o
menos este. Es el primer código de R medianamente largo que hago así que
habrá muchas cosas a mejorar o corregir así que no lo tomes ni mucho menos
como una referencia.
#Inicializamos variables
AnoModelo <- 2010
Tipologia <- 1
#Cargamos la libreria RODBC y forecast para calcular modelos ARIMA
automaticos
library(RODBC)
library(forecast)
#Establecemos la conexion con el SQL Server
channel <- odbcDriverConnect( "case=nochange; Driver=SQL Server;
Server=+++++++; Database=+++++++++; uid=++++++++; pwd=++++++++++;
wsid=++++++++++++;")
#Una vez establecida la conexion consultamos la tabla de la BBDD y la
guardamos en un data frame llamado Total
Total <- sqlQuery(channel,paste("SELECT * FROM
DATOS_FOMENTO_PROVINCIAL_BASE2005 WHERE codprovi is not null "))
#Definimos funcion para calcular p-valores de los coeficientes
cwp <- function (object)
{
coef <- coef(object)
if (length(coef) > 0) {
mask <- object$mask
sdev <- sqrt(diag(vcov(object)))
t.rat <- rep(NA, length(mask))
t.rat[mask] <- coef[mask]/sdev
pt <- 2 * pnorm(-abs(t.rat))
setmp <- rep(NA, length(mask))
setmp[mask] <- sdev
sum <- rbind(coef, setmp, t.rat, pt)
dimnames(sum) <- list(c("coef", "s.e.", "t
ratio", "p-valor"),
names(coef))
return(sum)
} else return(NA)
}
#Creamos df que contendra los datos del modelo
modelo <- data.frame(Tipologia = numeric(0),
Provincia = numeric(0),
AnoModelo = numeric(0),
ecuacion =character(0),
aicc = numeric(0),
bic =numeric(0),
sigma2 = numeric(0),
loglik = numeric(0),
coeficientes =
character(),
signific =character(),
MASE = numeric(0)
)
#Creamos df con datos de las predicciones
predicciones <- data.frame(Tipologia = numeric(0),
Provincia = numeric(0),
AnoModelo = numeric(0),
ecuacion =character(0),
Anio_Trimestre =
character(0),
estimacion = numeric(0),
inf8df.modelo0 =
numeric(0),
sup80 = numeric(0),
inf95 = numeric(0),
sup95 = numeric(0)
)
#Generamos bucle para recorrer las 52 provincias
for (i in 1:52)
{
Provincia <- i
#Generamos un df llamado activo donde guardaremos solo los datos
de la provincia que nos interesa del df Total
activo <- subset(Total, codprovi == i, select =
c(codprovi,Año,Trimestre,Precio_Totales))
#Generamos, a partir de activo, un objeto ts que llameremos st
para tener la serie temporal
st <- ts(activo$Precio_Totales, frequency = 4, start =
c(activo[1,2], activo[1,3]))
#Asignamos un objeto forecast llamado fit que contiene el mejor
modelo automaticamente ajustado a los datos
fit <- auto.arima(st, ic ="aicc", stepwise= TRUE)
#Obtenemos una variable con los coeficientes del modelo y sus
p-valores
pv <- as.data.frame(cwp(fit)) #Guardamos diagnosis de parametros
en un df
numparam <- length(colnames(pv)) #Obtenemos el numero de
parametros presente en el modelo
coeficientes <- data.frame( CF = character(0) ) #Inicializamos df
que contendra los coeficientes del modelo
signific <- TRUE #Hacemos el modelo valido por defecto
for (j in 1:numparam)
{
parame <- colnames(pv)[j] #Obtenemos los nombres de
los parametros que forman el modelo
coefic <- round(pv[1,j],4) #Obtenemos los
coeficientes de esos parametros
pvalor <- round(pv[4,j],4) #Obtenemos los p-valores
de los coeficientes (si p-valor < 0.1 --> parametro significativo)
if (parame != "cwp(fit)")
{
if ((pvalor >= 0.1) & ( signific = TRUE)
)
#Chequeamos los p-valores y si hay alguno > 0.1 marcamos
modelo como no valido
signific = FALSE
else
signific = TRUE
}
else
signific = FALSE
if (j == 1) #Construimos cadena donde se guardan los
nombres de los parametros, sus coeficientes y sus p-valores asociados
coeficientes <- paste("Param:" , parame,
"
Coefic:" , coefic, " P-Val: ", pvalor)
else
coeficientes <- paste(coeficientes, "|
Param:" ,
parame, " Coefic:" , coefic, " P-Val: ", pvalor)
}
#Extraemos los datos de fit que nos interesa guardar a un df para
añadirlos al df modelo
ecuacion <- paste("ARIMA("
,fit$arma[1],fit$arma[6],fit$arma[2],")(",fit$arma[3],
fit$arma[7],fit$arma[4], ")" )
MASE <-accuracy(fit)[6]
df.modelo <- data.frame(Tipologia, Provincia, AnoModelo, ecuacion,
fit$aicc , fit$bic, fit$sigma2, fit$loglik, coeficientes, signific, MASE)
names(df.modelo)[5] <- "aicc"
names(df.modelo)[6] <- "bic"
names(df.modelo)[7] <- "sigma2"
names(df.modelo)[8] <- "loglik"
#Acumulamos los datos del modelo en el df modelo
modelo <- rbind(modelo,df.modelo)
#Extraemos las predicciones a un df
df.predicciones <- data.frame(Tipologia, Provincia, AnoModelo,
ecuacion, row.names(data.frame(forecast(fit,
h=12))),data.frame(forecast(fit, h=12)))
names(df.predicciones)[5] <- "Anio_Trimestre"
names(df.predicciones)[6] <- "estimacion"
names(df.predicciones)[7] <- "inf80"
names(df.predicciones)[8] <- "sup80"
names(df.predicciones)[9] <- "inf95"
names(df.predicciones)[10] <- "sup95"
#Acumulamos los datos de las predicciones en el df predicciones
predicciones <- rbind(predicciones, df.predicciones)
#Imprimimos provincia por la que vamos en pantalla y refrescamos
pantalla
print (paste(i, ecuacion, signific))
flush.console()
#Generamos y guardamos plot con la serie
dev.new()
titulo <- paste("Modelo: ", ecuacion, " en la
provincia: ", i)
plot(forecast(fit,h=12), main = titulo )
savePlot(filename = paste(Tipologia, "-",i),type =
"jpg", device =
dev.cur(), restoreConsole = TRUE)
dev.off()
#Generamos y guardamos plot de diagnosis del modelo
dev.new()
tsdiag(fit)
savePlot(filename = paste("Diagnosis ",Tipologia,
"-",i),type =
"jpg", device = dev.cur(), restoreConsole = TRUE)
dev.off()
}
#Cambiamos los rownames del df modelo para que nos sirvan de clave
primaria de la tabla
rownames(modelo) <- paste(modelo$Tipologia, modelo$Provincia,
modelo$AnoModelo)
#Subimos el df modelo a la tabla FE_MODELO del servidor SQL
sqlSave(channel,modelo,tablename ="FE_MODELOS", addPK =TRUE, append =
TRUE
)
#Cambiamos los rownames del df predicciones para que nos sirvan de clave
primaria de la tabla
rownames(predicciones) <- paste(predicciones$Tipologia,
predicciones$Provincia, predicciones$AnoModelo,
predicciones$Anio_Trimestre)
#Subimos el df modelo a la tabla FE_PREDICCIONES del servidor SQL
sqlSave(channel,predicciones,tablename ="FE_PREDICCIONES", addPK
=TRUE,
append = TRUE )
#Cerramos la conexion
close(channel)
rm(channel)
Espero que te haya aclarado algunas cosas. Si necesitas algo mas no dudes
en comentarlo
Un saludo
From: justodejorge@gmail.com
To: Jluis GILSANZ/ES/EUROPE/GROUP@BNPPARIBAS
Date: 09/05/2011 11:44
Subject: Series temporales
Hola Jose Luis, he estado pendiente de los correos que has mantenido en la
lista. Curiosamente yo tengo un problema parecido al tuyo y me alegro que
iniciaras el hilo. Aunque llevo algunos años con R, soy nuevo en series.
Te rogaría si puedes ayudarme, en caso afirmativo si el código que has
generado podría ayudarme. Te comento mi situación.
Dirijo dos trabajos de investigación y en ambos estimamos la eficiencia
técnica de las empresas. En uno de ellos referente a líneas aéreas en el
período 2002-2009 con datos mensuales. El valor de la eficiencia técnica
oscila entre 0-1 y pretendemos una vez estimada la eficiencia aplicar un
modelo ARIMA. Hemos estimado con el comando auto.arima con y sin drift y
nos muestra estos resultados el gráfico que te adjunto.
Series: Eficiencia
ARIMA(1,1,1)
Call: auto.arima(x = Eficiencia, allowdrift = TRUE)
Coefficients:
ar1 ma1
0.7720 -0.9644
s.e. 0.0813 0.0363
sigma^2 estimated as 0.0001375: log likelihood = 287.15
AIC = -568.31 AICc = -568.04 BIC = -560.65
Antes te trabajar con R comence con statgraphics (es muy elemental, pero
tiene una ayuda que te indica como interpretar los resultados) y el modelo
que me indica después de contrastar es un ARMA (1,0). No acabo de
relacionar ambos resultados aunque entiendo que el válido es el
proporcionado con R aunque después de vuestros correos. Finalmente no se
si conoces algún trabajo publicado (los que encuentro son muy duros) que
sirva de referencia para conocer como explicar la estimación de una serie.
Espero no molestarte con este correo
Gracias de antemano
Saludos cordiales
Justo
--
***************************************************
Prof. Justo de Jorge Moreno.
Depart.de Ciencias Empresariales
Facultad de Ciencias Económicas y Empresariales.
Plaza de la Victoria 2
Alcalá de Henares 28802
Universidad de Alcalá
Teléfono: 918854293; Fax: 918854294
E-mail: justo.dejorge@uah.es
Web: www2.uah.es/justo_de_jorge
*************************************************[attachment
"serie.doc"
deleted by Jluis GILSANZ/ES/EUROPE/GROUP]
-- AVISO LEGAL --
Los datos personales que en esta comunicación aparecen, así como los que nuestra
empresa mantiene de Vd. y de su empresa, son tratados con la finalidad de
mantener
el contacto así como realizar las gestiones que en esta aparecen (Ley Orgánica
15/1999, de 13 de diciembre, de Protección de Datos de Carácter Personal).
Puede ejercer sus derechos de acceso, rectificación, cancelación y oposición
dirigiéndose a atencion.clientes@tasacionesh.com
La utilización de su dirección de correo electrónico por parte de nuestra
empresa
queda sujeta a las disposiciones de la Ley 34/2002, de Servicios de la Sociedad
de
la Información y el Comercio Electrónico. Si Vd. recibe comunicación comercial
por
nuestra parte y desea dejar de recibirla, rogamos nos lo comunique por vía
electrónica
a través de la dirección atencion.clientes@tasacionesh.com
[[alternative HTML version deleted]]
