Skip to Content

Fuzzy Hybrid System for Forecasting Time Series

pfv007's picture

Authors: Patricio Fuenmayor Viteri and Hermann Mena
Emails: patricio.fuenmayor@gmail.com
City: Quito - Ecuador

Introduction

The management's of scenarios in decision-making is a necessary tool for the treatment of uncertainty. These scenarios consider assumptions, estimates, models, etc.., wich are obtained based on technical and/or of business criteria. Technically, the choice of the best model is based on the optimization criterion chosen as reference information. This best model is used to forecasts in a given time horizon. In practice, the experts opinion do not necessarily agree with the estimated values. Hence the search for a model that combines expert opinion as well as statistically acceptable fitting becomes a necessary task. The idea of easing the model selection criterion provides a set of possible options which, in a global view, may show the tendency in the fitting and forecasting.

We propose to generate a list of linear models that are chosen based on search and discrimination criteria by minimizing the information criteria of Akaike and Schwarz. Based on this list we forecast using each model. To characterize the trend and pattern we use fuzzy clustering techniques which characterized the average forecast. The proposed methodology is a first step in the structuring of the so-called hybrid systems, which are characterized by complementation of two or more techniques together, in this case, automatic forecasting and fuzzy systems. Numerical experiments with the financial series of savings accounts of a financial institution are presented.

The source code of the system and the data are available at the url.

Fuzzy Hybrid System

The fuzzy hybrid system (FHS) consists of three main stages. The stages here built such that the data is structured. This strusture let therefore apply the algorithms which here implemented in specific routines. A plot funtion allows to visualize the results of each stage.

The first stage involves the generation of parameter matrices of the test models. The best model is chosen based on the minimization of information criteria (IC) using the auto.arima2 funtion. This is a variant of the auto.arima function in forecast package of R. Inorder to cover a wider range of models three criteria are used for information (Akaike or AIC, BIC and Bayesian or Akaike second order or AICc). The criteria is a imput of the function, then a parameter array is generated for each of them. A key aspect of this stage is the exclusion of the models that are considered by experts as ineligible forecasting models (since the behavior is not consistent with characteristics of the series, as are the trend, seasonality, etc.). These models are the result of the analysis and iterative process undertaken with the system. The models are stored and incorporated into the process as a filter, so that the resulting parameter list does not contain them.

In the second stage, forecasting with each model of the list is performed using the function forecast . The function requires as inputs: the data series, the orders of the parameters (p,d,q,P,D,Q) model, and if the model has a drift. The function generates data containing: the model estimates, adjustments, errors and forecasts. Also the results for the best model by AIC criterion. In addition, to these results list, sort structures for CI values are included. This are important elements in the next stage.

The final stage uses the forecasting matrix as a basis for analysis and inference. Forecasts can graphically show a pattern of behavior, so the fuzzy hybrid system is based on the identification of this pattern. Hence each month is analized separatelly, in order to find a value (mean forecast) representing the predictions of the model for this month. A pattern can be interpreted as the similarity or proximity of each element to each other. To find these sets we using the fuzzy clustering algorithm, also known as c-means , the algorithm allows a level of membership of elements to each set, making it ideal for issues of uncertainty. The algorithm requires some imput parameters, three centers are defined and therefore forecasts can be classified as: overestimation, estimation and underestimation. The generated sets and the number of elements from each let identify the pattern of each month.

The inference defined to obtain the median forecast is based on the implementation of two weight criteria for each cluster centers in the forecasts made for each month. The first criterion considers the number of elements in each grouping, the second criterion is based on the number of elements in each group from low IC models.

With lists of results for each stage, we proceed to make the graphs where results are displayed clearly understand: adjustments, errors, predictions, forecasts respective means and compared with the best model SARIMA.

As part of the validation of results, the test is performed with real data forecast three months, where the results of the series analyzed, validated the hypothesis for this work.

Results for saving account time series

For this section show the principal commands and results.

Source and libraries

source("newarima2_1.R") ## modified auto.arima funtion forecast package
source("shd.R")         ## hybrid fuzzy system
source("shd_aux.R")     ## auxiliar funtions (save graphs in png format)
source("utl.R")         ## utilities functions
 
require(forecast)
require(sqldf)
require(e1071)
require(RColorBrewer)
require(cwhmisc)
require(dummies)
require(corrplot)
require(gdata)
require(qualV)
require(xtable)

Read and setting data

# data
mes00 <- read.csv("mespch00.txt", header= TRUE, sep= "\t") 
# list of models to exclude
mdlexcaho <- read.csv("mdlexcaho.txt", header= TRUE, sep= "\t", colClasses=c("integer","character"))
 
# dimentions of data
mes_r <- dim(mes00)[1]  #number of rows
mes_c <- dim(mes00)[2]  #number of colunms
 
# define lenght of time series
mes_rp <- mes_r - 3
mes_rt <- mes_rp + 1
 
# test of forecastings
tstaho <- mes00[mes_rt:mes_r,4]

Funtions to hybrid fuzzy system

# define time series
dataho <- ts(mes00[1:mes_rp,4],start=c(2004,1),frequency=12)
# generate list of models
lstmdlaho <- shd_mdl(dataho,mdlexcaho,quaprb=1,trace=T)
# generate forecastings
lstprnaho <- shd_prn(lstmdlaho,hrzprn=12,nivitv=63,fil="rstprnaho.Rdata")
# fuzzy clustering and inference in the groups (mean forecast)
lstcdpaho <- shd_cdp(lstprnaho,fil="rstcdpaho.Rdata")
# generate graphs and save in pdf format
shd_grf(lstprnaho,lstcdpaho,tstaho,pre="aho",uni="millones de ")
# generate table results
shd_tbl(lstprnaho,lstcdpaho,tstaho,pre="aho")

Graphs and Results

The saving accounts time series and ACF, PACF functions
This help to analize the time serie,

List of models testeds

> lstmdlaho$mtzmdl  # list of models
   p d q P D Q c     mdl      aic      bic     aicc
1  2 0 2 1 1 1 1 2021111 814.9766 833.9322 817.0338
2  1 0 0 1 1 0 1 1001101 813.0464 822.5242 813.5869
3  1 0 0 0 1 0 1 1000101 818.1071 825.2154 818.4271
4  1 0 0 2 1 0 1 1002101 813.5041 825.3514 814.3261
5  1 0 0 1 1 1 1 1001111 813.8907 825.7379 814.7126
6  2 0 0 1 1 0 1 2001101 810.8385 822.6857 811.6604
7  2 0 1 1 1 0 1 2011101 812.4894 826.7061 813.6561
8  3 0 1 1 1 0 1 3011101 812.3404       NA 813.9179
9  2 0 0 0 1 0 1 2000101 817.2048       NA 817.7454
10 2 0 0 2 1 0 1 2002101 811.0993       NA 812.2660
11 2 0 0 1 1 1 1 2001111 811.7167       NA 812.8833
12 3 0 0 1 1 0 1 3001101 812.0306       NA 813.1973
13 1 0 1 1 1 0 1 1011101       NA 822.5266       NA

The adjustments of models (tendency)

The residuals and relative residuals of models
This show that the adjustments of models have a tendency and a pattern.

The residuals boxplot in the time. This show the inusual data and the limitations of the models to explain the time series

Correlation of forcastings

The forecastings of models and inference (mean forecast)
The first graph show the forcastings od models (the green line is the best SARIMA model).
The second graph show the inference (black line) The others models are draws in blue color, the intensity is major when the IC of model is lower.

Test of model with real data and confidense intervals
This show the shd model (black), the best ARIMA (green), the real values (red).

Conclusions and Recommendations

The work shows an encouraging perspective on the use of computational intelligence in dealing with financial data. In particular, the principles of fuzzy logic as a basis for the development of clustering techniques and their incorporation into the hybrid system, as they allow financial phenomena modeling under uncertainty. It should be noted that the methodology is not restricted to applications in finance, but can be used in various fields.

The modeling technique based on fuzzy hybrid systems can be understood as a local method, which performs a fuzzy partition of the domain into a finite number of regions, in which inferences are applied to identify patterns. These patterns allowing those clusters and obtain a value representative of the entire domain, thus unable to understand the behavior of appropriate models in a joint.

The proposed method is computationally efficient regarding processing time and memory requirements, this allows large-scale systems modeling (hundreds of sets with a history of more than 10 years with monthly data).

The use of forecasts allows automatic adjustments to find models that provide statistically acceptable. The development and model-specific phenomenon, with more complex techniques, for which the results can be considered as the starting point of this analysis, is an open door for each specific application.

In the field of business, it is important to use all kinds of scenarios that give an estimate of how you can involve a specific analysis on this basis and make the right decision. Thus, you should not ignore the view of experts, because their experience is valuable for understanding the phenomena, and puden become deciding factors in conditions of uncertainty. The inclusion criteria of experts is a great advantage of the proposed method.

The use of R as a deployment tool provides a great advantage in terms of: data management, visualization of results and implementation of methodologies, and functions that disponde incorporating specific analysis and have been used and validated by the scientific community uses this tool.

The use of scripts leaves open the opportunity to further develop the technique and incorporating more criteria, for example, the incorporation of modeling the variance, or even the ability to define new criteria of inference, which can give the shift to more robust estimates.

References

We are working in the english version of the paper, Fuzzy Hybrid System for Forecasting Time Series by Patricio Fuenmayor Viteri and Hermann Mena (preprint).

If you are interested in this work, please write to the e-mail.

Annexes (part of source code)

# FUNCION: shd_mdl
#   genera lista de objetos resultantes de la modelación automática: matriz de modelos y mejores  modelos en base a los criterios de información (CI)
# PARAMETROS:
#   dat: datos o serie financiera
#   mdlexc: modelos excluídos (archivo dat\mdlexcxxx.txt)
#   quaprb: probabilidad para seleccionar los modelos con ic menor al quantile(probs=quaprb) 0=min, 1=max (default)
#   trace: para imprimir en consola los modelos testeados. Por defecto FALSE
#   fil: nombre del archivo para guardar resultados (por defecto rstmdl.Rdata)
# SALIDA (lstmdl):
#   mtzmdl: matriz de modelos
#   mdl_aic: mejor modelo en base a aic
#   mdl_bic: mejor modelo en base a bic
#   mdl_aicc: mejor modelo en base a aicc
 
shd_mdl <- function(dat,mdlexc,quaprb=1,trace=FALSE,fil="rstmdl.Rdata")
{
setwd(dirwrk)
 
# GENERAR MODELOS
mdl_aic <- auto.arima2(dat,ic="aic",trace=trace,file="rst_aic.Rdata")
mdl_bic <- auto.arima2(dat,ic="bic",trace=trace,file="rst_bic.Rdata")
mdl_aicc <- auto.arima2(dat,ic="aicc",trace=trace,file="rst_aicc.Rdata")
 
# CARGAR RESULTADOS
load("rst_aic.Rdata")
rst_aic <- with(rst,data.frame(rst,icn="aic",mdl=paste(p,d,q,P,D,Q,c,sep="")))
rst_aic <- rst_aic[rst_aic$ic<1e20,]
 
load("rst_bic.Rdata")
rst_bic <- with(rst,data.frame(rst,icn="bic",mdl=paste(p,d,q,P,D,Q,c,sep="")))
rst_bic <- rst_bic[rst_bic$ic<1e20,]
 
load("rst_aicc.Rdata")
rst_aicc <- with(rst,data.frame(rst,icn="aicc",mdl=paste(p,d,q,P,D,Q,c,sep="")))
rst_aicc <- rst_aicc[rst_aicc$ic<1e20,]
 
rst_aic <- rst_aic[rst_aic$ic <= quantile(rst_aic$ic,probs=quaprb),]
rst_bic <- rst_bic[rst_bic$ic <= quantile(rst_bic$ic,probs=quaprb),]
rst_aicc <- rst_aicc[rst_aicc$ic <= quantile(rst_aicc$ic,probs=quaprb),]
 
mtzmdl <- droplevels(rbind(rst_aic,rst_bic,rst_aicc))
mtzmdl <- sqldf("select * from mtzmdl where mdl not in (select t01.mdl from mdlexc t01)")
 
mtzmdl01 <- mtzmdl[,c(-8,-9)]
mtzmdl01 <- sqldf("select distinct * from mtzmdl01")
 
mtzmdl02 <- data.frame(aggregate.table(mtzmdl$ic, mtzmdl$mdl, mtzmdl$icn, sum))
mtzmdl02 <- data.frame('mdl'=rownames(mtzmdl02),mtzmdl02)
 
mtzmdl <- sqldf("select t01.*, t02.aic, t02.bic, t02.aicc from mtzmdl01 t01 left join mtzmdl02 t02 on t01.mdl = t02.mdl")
names(mtzmdl) <- c("p","d","q","P","D","Q","c","mdl","aic","bic","aicc")
 
# GUARDAR MATRIZ DE MODELOS
save(mtzmdl,file=fil)
 
return(list("mtzmdl"=mtzmdl,"mdl_aic"=mdl_aic,"mdl_bic"=mdl_bic,"mdl_aicc"=mdl_aicc))
}
 
#------------------------------------------------------------------------------------------------------------
# FUNCION: shd_prn
#   genera lista con matrices de: ajustes, residuos y pronósticos. A más de esto, guarda los resultados del mejor modelo obtenido en base al AIC
# PARAMETROS:
#   lstmdl: lista de modelación
#   hrzprn: horizonte de pronóstico (h en la función forecast)
#   nivitv: nivel de confianza de los intervalos en la función forecast
#   nivsel: nivel se selección de los modelos (valor del aic max)
#   fil: nombre del archivo para guardar resultados. Por defecto (rstprn.Rdata)
# SALIDA (lstprn)
#   mtzaju: matriz de ajustes
#   mtzrsd: matriz de residuos
#   mtzrel: matriz de residuos relativos
#   mtzprn: matriz de pronósticos
#   mtzaic: matriz de valores de AIC
#   mtzajumej: matriz de ajuste del mejor modelo
#   mtzrsdmej: matriz de residuos del mejor modelo
#   mtzrelmej: matriz de residuos relativos del mejor modelo
#   mtzprnmej: matriz de pronósticos del mejor modelo
#   dat: datos o serie financiera
#   mtzmdl: matriz de modelos
 
shd_prn <- function(lstmdl,hrzprn=12,nivitv=63,nivsel=1e20,fil="rstprn.Rdata")
{
setwd(dirwrk)
 
mtzmdl <- lstmdl$mtzmdl
dat <- lstmdl$mdl_aic$x
 
# DIMENSIONES DE LAS MATRICES
 
nummdl <- dim(mtzmdl)[1]
numdat <- length(dat)
prd <- frequency(dat)
 
lstprn <- list(NULL)
 
mtzaju <- matrix(NA,nrow=numdat,ncol=nummdl)  # matriz ajustes
mtzrsd <- mtzaju                              # matriz residuos
mtzrel <- mtzaju                              # matriz residuos relativos
mtzprn <- matrix(NA,nrow=hrzprn,ncol=nummdl)  # matriz pronósticos
mtzaic <- matrix(NA,nrow=1,ncol=nummdl)       # matriz criterio de información AIC
 
for(i in 1:nummdl){
# CARGAR PARAMETROS DEL MODELO
mdlpar <- mtzmdl[i,]
mdl <- with(mdlpar,Arima(dat,order=c(p,d,q),seasonal=list(order=c(P,D,Q),period=prd),include.drift= c==1))
lstprn[[i]] <- forecast(mdl,h=hrzprn,level=nivitv)
if(lstprn[[i]]$model$aic <= nivsel)
{
mtzaju[,i] <- lstprn[[i]]$fit
mtzrsd[,i] <- lstprn[[i]]$residuals
mtzrel[,i] <- lstprn[[i]]$residuals/lstprn[[i]]$x
mtzprn[,i] <- lstprn[[i]]$mean
mtzaic[,i] <- lstprn[[i]]$model$aic
}
}
# ELIMINA COLUMNAS CON VALORES NA
mtzaju <- clean.na(mtzaju,margin=2)
mtzrsd <- clean.na(mtzrsd,margin=2)
mtzrel <- clean.na(mtzrel,margin=2)
mtzprn <- clean.na(mtzprn,margin=2)
mtzaic <- clean.na(mtzaic,margin=2)
 
# IDENTIFICACIÓN POR AIC
col01 <- sort(brewer.pal(3,"Blues"),decreasing=F) # display.brewer.pal(3, "Blues")
 
minaic <- min(mtzaic)
maxaic <- max(mtzaic)
itvaic <- (maxaic-minaic)/3
 
mtzaic<- t(mtzaic)
mtzaic <- data.frame("aic"=mtzaic,"col"=cut(mtzaic,include.lowest=T,breaks=seq(from=minaic,to=maxaic,by=itvaic), right=F,labels=c(col01)))
 
# CARGAR PARAMETROS DEL MEJOR MODELO SEGUN AIC
prnmej <- forecast(lstmdl$mdl_aic,h=hrzprn,level=nivitv)
mtzajumej <- prnmej$fit
mtzrsdmej <- prnmej$residuals
mtzrelmej <- prnmej$residuals/prnmej$x
mtzprnmej <- prnmej$mean
 
# GUARDAR LISTA DE PRONOSTICOS
save(lstprn,file=fil)
 
return(list("mtzaju"=mtzaju,"mtzrsd"=mtzrsd,"mtzrel"=mtzrel,"mtzprn"=mtzprn,"mtzaic"=mtzaic,"prnmej"=prnmej,"mtzajumej"=mtzajumej,"mtzrsdmej"=mtzrsdmej,"mtzrelmej"=mtzrelmej,"mtzprnmej"=mtzprnmej,"dat"=dat,"mtzmdl"=mtzmdl))
}
 
#------------------------------------------------------------------------------------------------------------
# FUNCION: shd_cdp
#   genera agrupamientos difusos de los pronósticos
# PARAMETROS:
#   lstprn: lista de ajustes y pronósticos
#   tipifr: tipo de inferencia (max = maximo, ptc = participación)
#   fil: nombre del archivo para guardar resultados. Por defecto (rstcdp.Rdata)
# SALIDA (lstcdp):
#   lstclu: lista de resultados del algoritmo c-medias
#   mtzclu: matriz de clusters de los pronósticos
#   mtzcen: matriz de centros
#   mtzele: matriz de número de elementos por grupo
#   mtzifr: matriz resultante de la inferencia
 
shd_cdp <- function(lstprn,tipifr="ptc",fil="rstcdp.Rdata"){
setwd(dirwrk)
 
mtzprn <- lstprn$mtzprn
numprn <- dim(mtzprn)[1]
nummdl <- dim(mtzprn)[2]
lstclu <- list(NULL)
 
mtzclu <- matrix(NA,nrow=nummdl,ncol=numprn)  # matriz clusterns
mtzcen <- matrix(NA,nrow=3,ncol=numprn)       # matriz centros
mtzele <- mtzcen                              # matriz elementos
mtzifr01 <- matrix(NA,nrow=1,ncol=numprn)     # matriz inferencia
 
for(i in 1:numprn){
# CREAR CLUSTERS PARA CADA VALOR DE PRONOSTICO, SI NO SE PUEDE REALIZAR EL CLUSTER, GRABA EL VALOR PROMEDIO
if (sd(mtzprn[i,]) <= 5)
lstclu[[i]] <- list("cluster"=rep(c(1,2,3),length.out=nummdl),"centers"=rep(mean(mtzprn[i,]),3))
else
lstclu[[i]] <- cmeans(mtzprn[i,],3,1000,verbose=F,method="cmeans")
 
mtzclu[,i] <- lstclu[[i]]$cluster
mtzcen[,i] <- lstclu[[i]]$centers
mtzele[,i] <- table(mtzclu[,i])
}
 
# TIPOS DE INFERENCIA (máximo o participación)
# CRITERIO DE PARTICIPACIÓN
if(tipifr=="max")
{
posmax <- max.col(t(mtzele))
for(i in 1:numprn){
mtzifr01[i] <- mtzcen[posmax[i],i]
}
}
else 
{
mtzptc01 <- prop.table(mtzele,margin=2)
mtzifr01 <- apply(mtzcen*mtzptc01,2,sum)
}
 
# PONDERA POR EL NIVEL DE AJUSTE DEL MODELO
dumcol <- data.frame(dummy(lstprn$mtzaic$col))
numcol <- dim(dumcol)[2]                        #variable reutilizada en el for
names(dumcol) <- paste("C",1:numcol,sep="")
mtzcol <- data.frame(mtzclu)
mtzcol <- data.frame(mtzcol,dumcol)
 
mtzacu <- matrix(NA,nrow=3,ncol=numprn)
nomvar <- as.vector(names(mtzcol))
 
for(i in 1:numprn){
var <- nomvar[i]
numcol <- fn$sqldf("select $var, sum(C1) as 'C' from mtzcol group by $var")
mtzacu[,i] <- numcol[,2]
}
 
mtzptc02 <- prop.table(mtzacu,margin=2)
mtzifr02 <- apply(mtzcen*mtzptc02,2,sum)
 
mtzifr <- rbind(mtzifr01,mtzifr02)
mtzifr <- apply(mtzifr,2,mean)
 
# GUARDAR RESULTADOS DE CLUSTERS DIFUSOS DE PRONOSTICOS
 
rstcdp <- list("lstclu"=lstclu,"mtzclu"=mtzclu,"mtzcen"=mtzcen,"mtzele"=mtzele,"mtzifr"=mtzifr)
save(rstcdp,file=fil)
 
return(rstcdp)
}
 
#------------------------------------------------------------------------------------------------------------
# FUNCION: shd_grf
#   genera gráficas 
# PARAMETROS:
#   lstprn: salida de shd_prn
#   lstcdp: salida de shd_cdp
#   tst: datos para test
#   pre: prefijo para nombre de gráficas
#   uni: unidades (leyendas)
# SALIDAS:
#   _dat.pdf: gráfica de la serie de datos y de las funciones ACF y PACF
#   _datdif.pdf: gráfica de la serie de datos diferenciada
#   _lag.pdf: gráfica de los retardos de la serie
#   _prn.pdf: gráfica de los pronósticos del mejor modelo
#   _rsd.pdf: gráfica de los residuos del mejor modelo
#   _dcpstl.pdf: gráfica de descomposición stl (método loess)
#   _dcpmul.pdf: gráfica de descomposición multiplicativa
#   _dcpadd.pdf: gráfica de descomposición aditiva
#   _dcpmul_s.pdf: gráfica del componente estacional multiplicativo
#   _dcpadd_s.pdf: gráfica del componente estacional aditivo
#   _dcpadd_sr.pdf: gráfica del componente estacional y del residuo (aditivos)
#   _mtzaju.pdf: gráfica de la matriz de ajustes
#   _mtzrsd.pdf: gráfica de la matriz de residuos
#   _mtzifr.pdf: gráfica de la matriz de pronósticos e inferencia
#   _mtztst.pdf: gráfica de la matriz de inferencia con pruebas de pronósticos e identificadores de los modelos
#   _mtzdat.pdf: gráfica de la matriz de datos conjuntos (dat + prn + inf)
#   _boxrsd.pdf: gráfica en boxplot de los residuos
#   _boxrel.pdf: gráfica en boxplot de los residuos relativos
#   _disrsd.pdf: gráfica de dispersión de los residuos y residuos relativos por modelo
 
shd_grf <- function(lstprn,lstcdp,tst,pre="",uni="")
{
mtzdat <- rbind(lstprn$mtzaju,lstprn$mtzprn)
mtzdatmej <- rbind(as.matrix(lstprn$mtzajumej),as.matrix(lstprn$mtzprnmej))
 
numdat <- length(lstprn$dat)
numprn <- numdat + 1
hrzprn <- dim(lstprn$mtzprn)[1]
numtot <- length(mtzdatmej)
col01 <- as.character(lstprn$mtzaic$col)
nummdl <- dim(lstprn$mtzprn)[2]
ordmdl <- 1:nummdl
 
low <- lstcdp$mtzifr*(1-(0.03*sqrt(1:hrzprn)))
hgh <- lstcdp$mtzifr*(1+(0.03*sqrt(1:hrzprn)))
 
labx <- rep(hrzprn,nummdl)
laby <- lstprn$mtzprn[hrzprn,]
miny <- min(min(lstprn$mtzprn),min(tst),min(low)) - 5
maxy <- max(max(lstprn$mtzprn),max(tst),max(hgh)) + 5
 
# DESCOMPOSICIÓN
dcpadd <- decompose(lstprn$dat,type="additive")
dcpmul <- decompose(lstprn$dat,type="multiplicative")
 
# GRAFICACIÓN
setwd(dirimg)
win.graph(width=11,height=7,pointsize=10)
 
tsdisplay(lstprn$dat,points=F,main="",xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
dev.copy2pdf(file=paste(pre,"_dat.pdf",sep=""))
 
tsdisplay(diff(lstprn$dat),points=F,main="",xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
dev.copy2pdf(file=paste(pre,"_datdif.pdf",sep=""))
 
lag.plot(lstprn$dat,4)
dev.copy2pdf(file=paste(pre,"_lag.pdf",sep=""))
 
plot(lstprn$prnmej,ylab=paste("saldos (",uni,"dólares)",sep=""))
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
dev.copy2pdf(file=paste(pre,"_prn.pdf",sep=""))
 
tsdisplay(lstprn$prnmej$residuals,points=F,main="",xlab="tiempo (meses)",ylab=paste("residuos (",uni,"dólares)",sep=""))
dev.copy2pdf(file=paste(pre,"_rsd.pdf",sep=""))
 
plot(stl(lstprn$dat,s.window=12))
dev.copy2pdf(file=paste(pre,"_dcpstl.pdf",sep=""))
 
plot(dcpmul)
dev.copy2pdf(file=paste(pre,"_dcpmul.pdf",sep=""))
 
plot(dcpadd)
dev.copy2pdf(file=paste(pre,"_dcpadd.pdf",sep=""))
 
plot(dcpmul$figure,type="l",xlab="tiempo (meses)",ylab="estacionalidad (factores)")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
dev.copy2pdf(file=paste(pre,"_dcpmul_s.pdf",sep=""))
 
plot(dcpadd$figure,type="l",xlab="tiempo (meses)",ylab=paste("estacionalidad (",uni,"dólares)",sep=""))
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
dev.copy2pdf(file=paste(pre,"_dcpadd_s.pdf",sep=""))
 
plot(dcpadd$random,type="l",col="blue",xlab="tiempo (meses)",ylab=paste("estacionalidad y residuos (",uni,"dólares)",sep=""))
lines(dcpadd$seasonal,col="red")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("estacionalidad","residuos"),fill=c("red","blue"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_dcpadd_sr.pdf",sep=""))
 
plot(1:numdat,lstprn$dat,type="l",lty=1,col="red",xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
matpoints(lstprn$mtzaju,type="l",lty=1,col=col01)
points(1:numdat,lstprn$mtzajumej,type="l",lty=1,col="green")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("real","mejor","otros"),fill=c("red","green","steelblue2"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_mtzaju.pdf",sep=""))
 
corrplot(cor(lstprn$mtzprn), method="ellipse")
dev.copy2pdf(file=paste(pre,"_mtzcor.pdf",sep=""))
 
#-----------------------------------------------------------------------------------------------------------------
win.graph(width=11, height=7,pointsize=10)
layout(matrix(c(1,2), byrow=T, ncol=1))
matplot(lstprn$mtzrsd,type="l",lty=1,col=col01,main="Residuos",xlab="tiempo (meses)",ylab=paste(uni,"dólares",sep=""))
points(1:numdat,lstprn$mtzrsdmej,type="l",lty=1,col="green")
abline(a=0,b=0,col="black")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("mejor","otros"),fill=c("green","steelblue2"),xjust=0.15,yjust=1,cex=0.8)
 
matplot(lstprn$mtzrel*100,type="l",lty=1,col=col01,main="Residuos relativos", xlab="tiempo (meses)",ylab="porcentaje")
points(1:numdat,lstprn$mtzrelmej*100,type="l",lty=1,col="green")
abline(a=0,b=0,col="black")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("mejor","otros"),fill=c("green","steelblue2"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_mtzrsd.pdf",sep=""))
graphics.off()
 
#-----------------------------------------------------------------------------------------------------------------
win.graph(width=11, height=7,pointsize=10)
layout(matrix(c(1, 2), byrow=F, ncol=2))
matplot(lstprn$mtzprn,type="l",lty=1,col=col01,xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
points(1:hrzprn,lstprn$mtzprnmej,type="l",lty=1,lwd=2,col="green")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("mejor","otros"),fill=c("green","steelblue2"),xjust=0.15,yjust=1,cex=0.8)
 
matplot(lstprn$mtzprn,type="l",lty=1,col=col01,xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
points(1:hrzprn,lstcdp$mtzifr,type="l",lty=1,lwd=2,col="black")
points(1:hrzprn,lstprn$mtzprnmej,type="l",lty=1,lwd=2,col="green")
matpoints(t(lstcdp$mtzcen),type="p",lty=1,col="red",cex=0.8)
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("shd","mejor","otros"),fill=c("black","green","steelblue2"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_mtzifr.pdf",sep=""))
graphics.off()
 
#-----------------------------------------------------------------------------------------------------------------
win.graph(width=11, height=7,pointsize=10)
matplot(lstprn$mtzprn,type="n",ylim=c(miny,maxy),xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
polygon(c(1:hrzprn,rev(1:hrzprn)),c(low,rev(hgh)),col="grey87",border=F,fillOddEven=T)
matpoints(lstprn$mtzprn,type="l",lty=1,col=col01)
points(1:hrzprn,lstcdp$mtzifr,type="l",lty=1,lwd=2,col="black")
points(1:hrzprn,lstprn$mtzprnmej,type="l",lty=1,lwd=2,col="green")
points(tst,type="l",lty=1,lwd=2,col="red")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
text(labx,laby,ordmdl,cex=0.5)
legend("bottomright",legend=c("real","shd","mejor","otros","zona conf."),fill=c("red","black","green","steelblue2","grey87"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_mtztst.pdf",sep=""))
 
#-----------------------------------------------------------------------------------------------------------------
matplot(mtzdat,type="l",lty=1,col=col01,xlab="tiempo (meses)",ylab=paste("saldos (",uni,"dólares)",sep=""))
points(1:numtot,mtzdatmej,type="l",lty=1,col="green")
points(1:numdat,lstprn$dat,type="l",lty=1,col="red",xlab=NA)
points(numprn:numtot,lstcdp$mtzifr,type="l",lty=1,col="black",xlab=NA)
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
abline(v=numprn,col="black")
legend("bottomright",legend=c("real","shd","mejor","otros"),fill=c("red","black","green","steelblue2"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_mtzdat.pdf",sep=""))
 
#-----------------------------------------------------------------------------------------------------------------
boxplot(t(lstprn$mtzrsd),outpch=".",xlab="tiempo (meses)",ylab=paste("residuos (",uni,"dólares)",sep=""))
lines(1:numdat,lstprn$mtzrsdmej,type="l",lty=1,col="green")
lines(1:numdat,t(apply(lstprn$mtzrsd,1,median)),type="l",lty=1,col="orange")
abline(a=0,b=0,col="black")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("mejor","otros (mediana)"),fill=c("green","orange"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_boxrsd.pdf",sep=""))
#-----------------------------------------------------------------------------------------------------------------
 
boxplot(t(lstprn$mtzrel*100),outpch=".",xlab="tiempo (meses)",ylab="residuos (porcentaje)")
lines(1:numdat,lstprn$mtzrelmej*100,type="l",lty=1,col="green")
lines(1:numdat,t(apply(lstprn$mtzrel*100,1,median)),type="l",lty=1,col="orange")
abline(a=0,b=0,col="black")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
legend("bottomright",legend=c("mejor","otros (mediana)"),fill=c("green","orange"),xjust=0.15,yjust=1,cex=0.8)
dev.copy2pdf(file=paste(pre,"_boxrel.pdf",sep=""))
 
#-----------------------------------------------------------------------------------------------------------------
layout(matrix(c(1, 2), byrow=F, ncol=2))
par(las=1)
boxplot(lstprn$mtzrsd,horizontal = T,outpch=20,names=lstprn$mtzmdl$mdl,xlab=paste("residuos (",uni,"dólares)",sep=""),ylab=NA)
abline(a=0,b=0,col="black")
abline(v=0,col="black")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
 
boxplot(lstprn$mtzrel*100,horizontal = T,outpch=20,names=lstprn$mtzmdl$mdl,xlab="residuos relativos (porcentaje)",ylab=NA)
abline(a=0,b=0,col="black")
abline(v=0,col="black")
grid(nx=NULL,ny=NULL,col="lightgray",lty="dotted",lwd=par("lwd"),equilogs=TRUE)
dev.copy2pdf(file=paste(pre,"_disrsd.pdf",sep=""))
 
graphics.off()
setwd(dirwrk)
}
 
#------------------------------------------------------------------------------------------------------------
# FUNCION: shd_tbl
#   genera y guarda tablas
# PARAMETROS:
#   lstprn: salida de shd_prn
#   lstcdp: salida de shd_cdp
#   tst: datos para test
#   pre: prefijo para nombre de gráficas
# SALIDAS:
#   prnrea.txt: pronósticos y datos reales
#   t01.tex: tabla de estadísticos en latex
 
shd_tbl <- function(lstprn,lstcdp,tst,pre="")
{
setwd(dirwrk)
 
hrzprn <- dim(lstprn$mtzprn)[1]
hrztst <- length(tst)
tst <- matrix(data=tst,nrow=hrztst,ncol=1)
fil <- paste(dirwrk,"/",pre,"_prnrea.txt",sep="")
 
tst01 <- data.frame("ord"=1:hrzprn,"mej"=c(round(lstprn$mtzprnmej,0)),"shd"=round(lstcdp$mtzifr,0))
tst02 <- data.frame("ord"=1:hrztst,"rea"=tst)
 
res01 <- sqldf("select t01.*, t02.rea from tst01 t01 left join tst02 t02 on t01.ord=t02.ord")
res02 <- res01[1:hrztst,]
 
res03m <- data.frame("MDL"="mej","RMSE"=RMSE(res02[,4],res02[,2],type="dissimilarity"),
  "MAE"=MAE(res02[,4],res02[,2],type="dissimilarity"),
  "MAPE"=MAPE(res02[,4],res02[,2],type="dissimilarity"))
 
res03s <- data.frame("MDL"="shd","RMSE"=RMSE(res02[,4],res02[,3],type="dissimilarity"),
  "MAE"=MAE(res02[,4],res02[,3],type="dissimilarity"),
  "MAPE"=MAPE(res02[,4],res02[,3],type="dissimilarity"))
 
res03 <- rbind(res03m,res03s)
 
write.table(res01,file=fil,row.names=F,quote=F,sep="\t")
save.tex(xtable(res03,digits=2,align=c("c","c","r","r","r")),paste(dirwrk,"/",pre,"_t01.tex",sep=""))
}

Comments

lowell.it1's picture

Our flagship link building service – Paint It White. Links from manually created beautiful web 2.0 posts on premium blogs. Supported by an array of tier 2 and 3 links for maximum ranking increases
link building service