#con los datos del INE
ex<-as.matrix(read.csv2("H91_09_EXPO.csv", header=TRUE, sep=";", dec=","))
head(ex)
dx<-as.matrix(read.csv2("H91_09_DEAD.csv", header=TRUE, sep=";", dec=","))
#pesos
wx=ex*0+1
head(wx)
dim(wx)


qdata.ES<-structure(list(x=as.integer(c(1:101)), y=as.integer(c(1991:2009)),
etx=structure(c(t(ex)), .Dim = as.integer(c(19,101))), dtx=structure(c(t(dx)),
.Dim = as.integer(c(19,101))), wa=structure(c(t(wx)),  .Dim = as.integer(c(19,101))),
.Names = c("x", "y", "etx", "dtx", "wa")))


 x=qdata.ES$x
 y=qdata.ES$y
 etx=qdata.ES$etx
 dtx=qdata.ES$dtx
 wa=qdata.ES$wa



#Cargando las funciones
source("fitModels.r")
source("simModels.r")


#modelo M1
resM1=fit701(x,y,etx,dtx,wa)

#Modelo M5
resM5=fit705(x,y,etx,dtx,wa)

#Modelo M7
resM7=fit707(x,y,etx,dtx,wa)


#Bayesian information criterion
resM1$BIC
resM5$BIC
resM7$BIC


#Objetos de cada modelo
names(resM1)
names(resM5)
names(resM7)


#Resultados M1
#GRáfico 3D
persp(resM1$y,resM1$x,log(resM1$mhat), phi=30, theta=55,col = "green3")
#Gráfico de los parámetros estimados para el modelo M1
par(mfrow=c(3,2))
 plot(resM1$cy,resM1$gamma3)
 plot(resM1$x,resM1$beta1)
 plot(resM1$x,resM1$beta2)
 plot(resM1$x,resM1$beta3)
 plot(resM1$y,resM1$kappa2)



# Simulaciones
sresM1=sim2001(resM1$x,resM1$y,resM1$beta1, resM1$beta2,resM1$kappa2,nsim=1000,nyears=10,tmax=40,x0=60)

## sres$xv is a vector of ages taken by a cohort initially aged x0 in year sres$y[1]
#sres$tpa is a 2-dimensional array of simulated values for the cohort survival
#index S(t,x0). sres$tpa[i,j] is the simulated value of S(t,x0) in year y[i] and for
#sample path j out of nsim sample paths.
fan(sresM1$xv,sresM1$tpa)


# sres$y is a vector of future calendar years in the projection
#sres$qaa is a 3-dimensional array of simulated values for the random future
#mortality rates q(t,x). sres$qaa[i,j,k] is the simulated mortality rate for age
#sres$xx[i], calendar year sres$y[j] and for sample path k out of nsim sample
#paths.
x11()
fan(sresM1$y,sresM1$qaa[25,,])



#Resultados M5
resM5=fit705(x,y,etx,dtx,wa)
names(resM5)

#objetos
resM5$kappa1; resM5$kappa2; resM5$kappa4; resM5$beta2; reM5$beta3;
resM5$beta4; resM5$gamma3; resM5$x0; resM5$x; resM5$y;
resM5$cy; resM5$epsilon; resM5$mhat; resM5$mtx; res$M5napr; resM5$ll; resM5$BIC; resM5$wa

#Gráfico de los parámetros estimados para el modelo M1
 par(mfrow=c(3,2))
 plot(resM5$cy,resM5$gamma3)
 plot(resM5$y,resM5$kappa1)
 plot(resM5$x,resM5$beta2)
 plot(resM5$y,resM5$kappa2)
 plot(resM5$x,resM5$beta4)
 plot(resM5$y,resM5$kappa4)

sresM5=sim2005(resM5$x,resM5$y,resM5$kappa1, resM5$kappa2,nsim=1000,nyears=10,tmax=40,x0=60)

fan(sresM5$xv,sresM5$tpa)
x11()
fan(sresM5$y,sresM5$qaa[25,,])