#definir funciones necesarias
logit<-function(x){log(x/(1-x))}
epam<- function(x,y){mean(abs((x-y)/x))*100} #x observados, y ajustados
mse <- function(x,y){sqrt(mean((x-y)^2))}
zeta<- function(def,n,q){(def-n*q)/sqrt(n*q*(1-q))}
inv.logit<-function(x){exp(x)/(1+exp(x))}
logvero<-function(pob, def,q){def*log(q)+(pob-def)*log(1-q)}#logverosimilitud de una binomial, fijos pob y def se maximiza q
dev<-function(pob,def,p1,p2){2*sum(logvero(pob,def,p1))-2*sum(logvero(pob,def,p2))} #deviance


                    ###########################
                    #FUNCION DE LEE-CARTER#####
                    ###########################
#trans es la transformacion de las medidas de mortalidad: logit(qxt), log(muxt), log(mxt)
#pob la población para cada edad y año
#def las defunciones para cada edad y año
#m es la edad mínima de la que se parte
#t es el número de años para los que se realiza el ajuste
#n es el número de edades para las que se realiza el ajuste
#lo es 1 cuando utilizamos la transformacion logit y otro valor cuando usamos log
funcion.lc<-function(trans,m,t,n,pob,def,lo){
#los a_x son simplemente las medias a lo largo del tiempo de ln(m_{xt}), en general de la transformación de las medidas de mortalidad
                    ax<-tapply(trans,rep(seq(m,m+n-1),t),mean)
 
                    #podemos estimar, respectivamente, k_t y b_x como los primeros vectores derecho e izquierdo correspondientes al
                    #principal valor de la SVD de la matriz (trans-ax)
                    lee.carter<-svd(matrix(trans,nrow=n,ncol=t)-matrix(rep(ax,t),nrow=n,ncol=t), nu=n, nv=t)
                    kt<-lee.carter$d[1]*lee.carter$v[,1]
                    bx<-lee.carter$u[,1]
                    kt<-as.vector(sum(bx)*kt) #es muy importante que esté antes de transformar bx sino no multiplico y divido, sum(bx) por lo mismo
                    bx<-as.vector(bx/sum(bx))
                    if(kt[1] <0){kt<- -kt        #este paso se hace por si los kt crecen, entonces han de cambiarse los signos
                                bx<- -bx}
                    else{kt<- kt
                         bx<- bx }
                    #pero se ha obtenido ajustando los logaritmos no las probabilidades
                    #reestimacion de las kt
                    poblacion<-matrix(pob,nrow=n,ncol=t)
                    defunciones<-matrix(def,nrow=n,ncol=t)
                    
                    Dest<-vector("numeric", t)#estimaciones de las muertes segun el modelo
                    for(i in 1: t) Dest[i]<-sum((exp(ax+kt[i]*as.vector(bx))/(1+exp(ax+kt[i]*as.vector(bx))))*poblacion[,i])
                    
                    Dt<-vector("numeric", t)
                    for(i in 1: t) Dt[i]<-sum(defunciones[,i])#calcula el total de muertes para cada año
                    Pt<-vector("numeric", t)
                    for(i in 1: dim(poblacion)[2]) Pt[i]<-sum(poblacion[,i])#calcula el total de población para cada año
                    if(lo==1){
                    reestimar<-function(indice){ Dest2<-vector("numeric", t)
                    for(i in 1: t) Dest2[i]<-sum((exp(ax+indice[i]*as.vector(bx))/(1+exp(ax+indice[i]*as.vector(bx))))*poblacion[,i])
                        result<-sum((Dt-Dest2)^2)
                        result}}
                    else{
                    reestimar<-function(indice){ Dest2<-vector("numeric", t)
                    for(i in 1: t) Dest2[i]<-sum((exp(ax+indice[i]*as.vector(bx)))*poblacion[,i])
                        result<-sum((Dt-Dest2)^2)
                        result}}
                    lcmodelo<-nlm(reestimar,kt)
                    lcmodelo2<-nlm(reestimar,lcmodelo$estimate)
                    #VOY A RESTAR LA MEDIA DE KT PARA CONSEGUIR QUE SU SUMA SEA CERO.
                    kt<-lcmodelo2$est-mean(lcmodelo2$est)
                    ax<-ax+as.vector((bx*mean(lcmodelo2$est)))
                    #calculo de las estimaciones de las probabilidades, residuos y medidas de bondad
                    epam.year<-vector("numeric", t)  
                    mse.year<-vector("numeric", t) 
                    deviance.year<-vector("numeric", t)
                    lee     <-matrix(rep(ax,t), nrow=n, ncol= t)+(matrix(bx, nrow=n, ncol= 1)%*%matrix(kt, nrow=1, ncol= t))
                    res.svd  <-matrix(trans,nrow=n, ncol=t)-lee
                    if(lo==1){estimaciones<-inv.logit(lee)
                              for(i in 1: t){epam.year[i]<-epam(matrix(inv.logit(trans),nrow=n,ncol=t)[,i],estimaciones[,i])
                                             mse.year[i]<-mse(matrix(inv.logit(trans),nrow=n,ncol=t)[,i],estimaciones[,i])
                              deviance.year[i]<-dev(poblacion[,i],defunciones[,i],matrix(inv.logit(trans),nrow=n,ncol=t)[,i],estimaciones[,i])}}
                    else{estimaciones<-exp(lee)
                        for(i in 1: t){epam.year[i]<-epam(matrix(exp(trans),nrow=n,ncol=t)[,i],estimaciones[,i])
                                      mse.year[i]<-mse(matrix(inv.logit(trans),nrow=n,ncol=t)[,i],estimaciones[,i])}}
                    prueba<-list(ax=ax,bx=bx,kt=kt,lee=lee,est=estimaciones,resid=res.svd,epam.year=epam.year,mse.year=mse.year,deviance.year=deviance.year)
                    prueba
                    }
#leer datos                
datos.brutos<-read.table("c:/congresos_internacionales/ime2004_roma/datosq.txt",header=T)
summary(datos.brutos)
#seleccionar para edades y tiempo que se quiere realizar el ajuste, columnas de cada ajuste
hombres<-datos.brutos[(datos.brutos$edad>=0) & (datos.brutos$tiempo<=1999),c(1,2,5,6)] #hombres
summary(hombres)
attach(hombres)
hombrescv<-funcion.lc(log(homqx), 0,20,97,homriesgo,homriesgo*homqx,0)
#seleccionar para edades y tiempo que se quiere realizar el ajuste, columnas de cada ajuste
mujeres<-datos.brutos[(datos.brutos$edad>=0) & (datos.brutos$tiempo<=1999),c(1,2,9,10)]#mujeres
summary(mujeres)
attach(mujeres)

mujqx[mujqx==0]<-3.342805e-007 #(0.01/poblacion)
mujerescv<-funcion.lc(log(mujqx), 0,20,97,mujriesgo,mujriesgo*mujqx,0)


plot(seq(1980,1999,1),hombrescv$kt,type="n",xlab="age",ylab="kt", ylim=c(-20,20))
lines(seq(1980,1999,1),hombrescv$kt,col=4)
lines(seq(1980,1999,1),mujerescv$kt,col=2)
legend(1980,-10,c("Men","Women"),lty=1,col =2:3)



plot(seq(0,96,1),hombrescv$bx,type="n",ylab="bx",,xlab="age")
lines(seq(0,96,1),hombrescv$bx,col=4)
lines(seq(0,96,1),mujerescv$bx,col=2)
lines(seq(0,96,1),rep(0,97))
legend(60,0,c("Men","Women"),lty=1,col =c(4,2))

plot(seq(0,96,1),hombrescv$ax,type="n",ylab="ax",xlab="age",ylim=c(-9,-1))
lines(seq(0,96,1),hombrescv$ax,col=4)
lines(seq(0,96,1),mujerescv$ax,col=2)
legend(70,-8,c("Men","Women"),lty=1,col =c(4,2))