################ METODOS ESTADISTICOS BAYESIANOS #####################
################ PROFESOR: Luis E. Nieto Barajas #####################

######################### COMANDOS DE R ##############################

######################################################################
#Diagnóstico de convergencia

#Lectura de datos----------------------------------------------------

cadena<-read.table("c:/Temp/cadena.txt")

#Gráfica de diagnóstico----------------------------------------------
par(mfrow=c(2,2))
plot(cadena[,1],cadena[,2],type="l")
plot(cadena[,1],cumsum(cadena[,2])/cadena[,1],type="l")
acf(cadena[,2])
hist(cadena[,2])


######################################################################
#Ejercicio 3: Regresion lineal simple

#Lectura de datos----------------------------------------------------
calif<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/calificaciones.txt",header=T)

#Análisis exploratorio-----------------------------------------------
par(mfrow=c(1,1))
plot(calif[,2],calif[,1],xlab="MO",ylab="SP")

#Ajuste del modelo-----------------------------------------------
predic3<-read.table("c:/Temp/predic3.txt",header=T,row.names=1)
par(mfrow=c(1,1))
plot(calif[,1],predic3[,1],xlab="SP",ylab="SP^")
lines(calif[,1],calif[,1])
Syy<-(19)*var(calif[,1])
SCR<-(19)*var(predic3[,1])
R2<-SCR/Syy
R2

#Análisis de residuos-----------------------------------------------

res<-calif[,1]-predic3[,1]
par(mfrow=c(2,1))
plot(predic3[,1],res,xlab="SP",ylab="residuos")
plot(calif[,2],res,xlab="MO",ylab="residuos")

par(mfrow=c(2,2))
plot(res[1:19],res[2:20],xlab="residuos(-1)",ylab="residuos")
acf(res)
qqnorm(res)
boxplot(res)

rese<-(calif[,1]-predic3[,1])*sqrt(5.275)
par(mfrow=c(2,1))
plot(predic3[,1],rese,xlab="SP",ylab="residuos est")
plot(calif[,2],rese,xlab="MO",ylab="residuos est")

par(mfrow=c(2,2))
plot(rese[1:19],rese[2:20],xlab="residuos est(-1)",ylab="residuos est")
acf(rese)
qqnorm(rese)
boxplot(rese)


######################################################################
#Ejercicio 4: Regresion lineal múltiple

#Lectura de datos----------------------------------------------------
salarios<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/salarios.txt",header=T)

#Análisis exploratorio----------------------------------------------
par(mfrow=c(1,1))
pairs(salarios)
cor(salarios)
boxplot(salarios)

#Ajuste del modelo-----------------------------------------------
predic4<-read.table("c:/Temp/predic4.txt",header=T,row.names=1)
par(mfrow=c(1,1))
plot(salarios[,1],predic4[,1],xlab="Y",ylab="Y^")
lines(salarios[,1],salarios[,1])
Syy<-(23)*var(salarios[,1])
SCR<-(23)*var(predic4[,1])
R2<-SCR/Syy
R2

#Análisis de residuos-----------------------------------------------

res<-salarios[,1]-predic4[,1]
par(mfrow=c(2,2))
plot(predic4[,1],res,xlab="Y",ylab="residuos")
plot(salarios[,2],res,xlab="X1",ylab="residuos")
plot(salarios[,3],res,xlab="X2",ylab="residuos")
plot(salarios[,4],res,xlab="X3",ylab="residuos")

par(mfrow=c(2,2))
plot(res[1:23],res[2:24],xlab="residuos(-1)",ylab="residuos")
acf(res)
qqnorm(res)
boxplot(res)

rese<-(salarios[,1]-predic4[,1])*sqrt(0.3261)
par(mfrow=c(2,2))
plot(predic4[,1],rese,xlab="Y",ylab="residuos est")
plot(salarios[,2],rese,xlab="X1",ylab="residuos est")
plot(salarios[,3],rese,xlab="X2",ylab="residuos est")
plot(salarios[,4],rese,xlab="X3",ylab="residuos est")

par(mfrow=c(2,2))
plot(rese[1:23],rese[2:24],xlab="residuos est(-1)",ylab="residuos est")
acf(rese)
qqnorm(rese)
boxplot(rese)


######################################################################
#Ejercicio 5: Modelo lineal dinámico (Producción de leche)

#Lectura de datos----------------------------------------------------
milk<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/milk.txt",header=T)

#Análisis exploratorio-----------------------------------------------
par(mfrow=c(2,1))
plot(milk[,2],milk[,1],xlab="#Vacas",ylab="Produccion")
plot(1:13,milk[,1],xlab="Tiempo",ylab="Produccion")

#Ajuste del modelo estático (x[t])-----------------------------------
predic5<-read.table("c:/Temp/predic5.txt",header=T,row.names=1)
par(mfrow=c(2,1))
plot(milk[,1],predic5[,1],xlab="Y",ylab="Y.New")
lines(milk[,1],milk[,1])
plot(milk[,2],milk[,1],xlab="#Vacas",ylab="Produccion")
lines(milk[,2],predic5[,1])
lines(milk[,2],predic5[,4],lty=2)
lines(milk[,2],predic5[,6],lty=2)
Syy<-(12)*var(milk[,1])
SCR<-(12)*var(predic5[,1])
R2<-SCR/Syy
R2

#Ajuste del modelo estático (Tiempo)-----------------------------------
predic5<-read.table("c:/Temp/predic5.txt",header=T,row.names=1)
par(mfrow=c(2,1))
plot(milk[,1],predic5[,1],xlab="Y",ylab="Y.New")
lines(milk[,1],milk[,1])
plot(1:13,milk[,1],xlab="Tiempo",ylab="Produccion",ylim=c(110,140))
lines(1:13,predic5[,1])
lines(1:13,predic5[,4],lty=2)
lines(1:13,predic5[,6],lty=2)
Syy<-(12)*var(milk[,1])
SCR<-(12)*var(predic5[,1])
R2<-SCR/Syy
R2

#Ajuste del modelo estático (Tiempo + x[t])----------------------------
predic5<-read.table("c:/Temp/predic5.txt",header=T,row.names=1)
par(mfrow=c(2,1))
plot(milk[,1],predic5[,1],xlab="Y",ylab="Y.New")
lines(milk[,1],milk[,1])
plot(1:13,milk[,1],xlab="Tiempo",ylab="Produccion",ylim=c(110,140))
lines(1:13,predic5[,1])
lines(1:13,predic5[,4],lty=2)
lines(1:13,predic5[,6],lty=2)
Syy<-(12)*var(milk[,1])
SCR<-(12)*var(predic5[,1])
R2<-SCR/Syy
R2

#Ajuste del modelo dinámico lineal (x[t])------------------------------
predic5<-read.table("c:/Temp/predic5.txt",header=T,row.names=1)
par(mfrow=c(2,2))
plot(milk[-1,1],predic5[,1],xlab="Y",ylab="Y.New")
lines(milk[-1,1],milk[-1,1])
plot(milk[,2],milk[,1],xlab="#Vacas",ylab="Produccion",ylim=c(113,142))
lines(milk[-1,2],predic5[,1])
plot(1:13,milk[,1],xlab="Tiempo",ylab="Produccion",ylim=c(113,142))
lines(2:13,predic5[,1])
plot(1:13,milk[,1],xlab="Tiempo",ylab="Produccion",ylim=c(108,142))
lines(2:13,predic5[,5])
lines(2:13,predic5[,4],lty=2)
lines(2:13,predic5[,6],lty=2)

#Ajuste del modelo dinámico lineal (Sin covariables)-------------------
predic5<-read.table("c:/Temp/predic5.txt",header=T,row.names=1)
par(mfrow=c(2,1))
plot(milk[-1,1],predic5[,1],xlab="Y",ylab="Y.New")
lines(milk[-1,1],milk[-1,1])
plot(1:13,milk[,1],xlab="Tiempo",ylab="Produccion",ylim=c(94,150))
lines(2:13,predic5[,1])
lines(2:13,predic5[,4],lty=2)
lines(2:13,predic5[,6],lty=2)


######################################################################
#Ejercicio 6: Modelo lineal dinámico (Participación de mercado)

#Lectura de datos----------------------------------------------------
mercado<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/mercado.txt",header=T)

#Análisis exploratorio-----------------------------------------------
par(mfrow=c(2,2))
plot(1:104,mercado[,4],xlab="Tiempo",ylab="Share")
plot(mercado[,1],mercado[,4],xlab="CPROM",ylab="Share")
plot(mercado[,2],mercado[,4],xlab="PRICE",ylab="Share")
plot(mercado[,3],mercado[,4],xlab="OPROM",ylab="Share")

#(a) Ajuste del modelo estático--------------------------------------
predic6<-read.table("c:/Temp/predic6.txt",header=T,row.names=1)
par(mfrow=c(2,1))
plot(mercado[-1,4],predic6[,1],xlab="Y",ylab="Y.New")
lines(mercado[-1,4],mercado[-1,4])
plot(1:104,mercado[,4],xlab="Tiempo",ylab="SHARE",ylim=c(38.5,44.5))
lines(2:104,predic6[,1])
lines(2:104,predic6[,4],lty=2)
lines(2:104,predic6[,6],lty=2)

#(b) Ajuste del modelo dinámico--------------------------------------
predic6<-read.table("c:/Temp/predic6.txt",header=T,row.names=1)
par(mfrow=c(2,1))
plot(mercado[-1,4],predic6[,5],xlab="Y",ylab="Y.New")
lines(mercado[-1,4],mercado[-1,4])
plot(1:104,mercado[,4],xlab="Tiempo",ylab="SHARE",ylim=c(38.5,45))
lines(2:104,predic6[,5])
lines(2:104,predic6[,4],lty=2)
lines(2:104,predic6[,6],lty=2)


######################################################################
#Ejercicio 7: Tasas de mortalidad 

#Lectura de datos----------------------------------------------------
mortality<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/mortality.txt",header=T)

#Análisis exploratorio-----------------------------------------------
tme<-mortality[,2]/mortality[,3]
par(mfrow=c(2,2))
plot(mortality[,1],mortality[,2],xlab="T.Exposición",ylab="No.Muertes")
plot(mortality[,1],tme,xlab="T.Exposición",ylab="T.Mortalidad")
plot(mortality[,1],log(tme/(1-tme)),xlab="T.Exposición",ylab="logit(T.Mortalidad)")

#Ajuste del modelo --------------------------------------------------
predic7<-read.table("c:/Temp/predic7.txt",header=T,row.names=1)
devia7<-read.table("c:/Temp/devia7.txt",header=T,row.names=1)
par(mfrow=c(2,2))
plot(mortality[,1],mortality[,2],xlab="T.Exposición",ylab="No.Muertes")
points(mortality[,1],predic7[,1],pch=3)
plot(mortality[,1],tme,xlab="T.Exposición",ylab="T.Mortalidad")
points(mortality[,1],predic7[,1]/mortality[,3],pch=3)
plot(mortality[,1],devia7[,1],xlab="T.Exposición",ylab="Devianza")
abline(h=3.06/2/6)


######################################################################
#Ejercicio 8: Número de desastres 

#Lectura de datos----------------------------------------------------
desastres<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/desastres.txt",header=T)

#Análisis exploratorio-----------------------------------------------
par(mfrow=c(1,1))
plot(desastres[,1],desastres[,2],xlab="Año",ylab="No.Desastres",type="b")

# (a) Ajuste del modelo (tasa en función del tiempo)-----------------
tasas8<-read.table("c:/Temp/tasas8.txt",header=T,row.names=1)
par(mfrow=c(1,1))
plot(desastres[,1],desastres[,2],xlab="Año",ylab="No.Desastres",type="n")
lines(desastres[,1],desastres[,2],lty=2)
lines(desastres[,1],tasas8[,1])
lines(desastres[,1],tasas8[,4],lty=3)
lines(desastres[,1],tasas8[,6],lty=3)

# (b) Ajuste del modelo (tasa constante en dos períodos)-------------
tasas8<-read.table("c:/Temp/tasas8.txt",header=T,row.names=1)
par(mfrow=c(1,1))
plot(desastres[,1],desastres[,2],xlab="Año",ylab="No.Desastres",type="n")
lines(desastres[,1],desastres[,2],lty=2)
lines(desastres[,1],tasas8[,1])
lines(desastres[,1],tasas8[,4],lty=3)
lines(desastres[,1],tasas8[,6],lty=3)
#abline(v=1891)

# (b) Ajuste del modelo (modelo dinámico)-------------
tasas8<-read.table("c:/Temp/tasas8.txt",header=T,row.names=1)
par(mfrow=c(1,1))
plot(desastres[,1],desastres[,2],xlab="Año",ylab="No.Desastres",type="n")
lines(desastres[,1],desastres[,2],lty=2)
lines(desastres[,1],tasas8[,1])
lines(desastres[,1],tasas8[,4],lty=3)
lines(desastres[,1],tasas8[,6],lty=3)


######################################################################
#Ejercicio 9: Incidencia y tasas de mortalidad 

#Lectura de datos----------------------------------------------------
leucemia<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/leucemia.txt",header=T)

#Análisis exploratorio-----------------------------------------------
par(mfrow=c(2,1))
plot(leucemia[,4],xlab="Categoría",ylab="No.Muertes")
plot(leucemia[,4]/leucemia[,6],xlab="Categoría",ylab="Frec.Rel")

# Ajuste del modelo (independencia)----------------------------------
tasas91<-read.table("c:/Temp/tasas91.txt",header=T,row.names=1)
par(mfrow=c(1,1))
ymin=min(tasas91[,4])
ymax= max(tasas91[,6])
plot(leucemia[,4]/leucemia[,6]*100000,xlab="Categoría",ylab="Tasas",ylim=c(ymin,ymax))
points(tasas91[,1],pch=3)
for (i  in 1:8){segments(i,tasas91[i,4],i,tasas91[i,6])}

# Ajuste del modelo (independencia + jerárquico)-------------------
tasas92<-read.table("c:/Temp/tasas92.txt",header=T,row.names=1)
par(mfrow=c(2,1))
ymin=min(tasas91[,4])
ymax= max(tasas91[,6])
plot(leucemia[,4]/leucemia[,6]*100000,xlab="Categoría",ylab="Tasas",ylim=c(ymin,ymax))
points(tasas91[,1],pch=3)
for (i  in 1:8){segments(i,tasas91[i,4],i,tasas91[i,6])}
ymin=min(tasas91[,4])
ymax= max(tasas91[,6])
plot(leucemia[,4]/leucemia[,6]*100000,xlab="Categoría",ylab="Tasas",ylim=c(ymin,ymax))
points(tasas92[,1],pch=4)
for (i  in 1:8){segments(i,tasas92[i,4],i,tasas92[i,6])}



######################################################################
#Ejercicio 10: Probabilidades de reclamación 

#Lectura de datos----------------------------------------------------
reclama<-read.table("http://allman.rhon.itam.mx/~lnieto/index_archivos/reclama.txt",header=T)

#Análisis exploratorio-----------------------------------------------
par(mfrow=c(2,1))
plot(reclama[,1],xlab="Categoría",ylab="No.Reclamaciones")
plot(reclama[,1]/reclama[,2],xlab="Categoría",ylab="Frec.Rel")

# Ajuste del modelo (independencia)----------------------------------
predic101<-read.table("c:/Temp/predic101.txt",header=T,row.names=1)
par(mfrow=c(1,1))
ymin=min(predic101[,4])
ymax= max(predic101[,6])
plot(reclama[,1]/reclama[,2],xlab="Categoría",ylab="Probabilidades",ylim=c(ymin,ymax))
points(predic101[,1],pch=3)
for (i  in 1:10){segments(i,predic101[i,4],i,predic101[i,6])}

# Ajuste del modelo (independencia + jerárquico)-------------------
predic102<-read.table("c:/Temp/predic102.txt",header=T,row.names=1)
par(mfrow=c(2,1))
ymin=min(predic101[,4])
ymax= max(predic101[,6])
plot(reclama[,1]/reclama[,2],xlab="Categoría",ylab="Probabilidades",ylim=c(ymin,ymax))
points(predic101[,1],pch=3)
for (i  in 1:10){segments(i,predic101[i,4],i,predic101[i,6])}
ymin=min(predic101[,4])
ymax= max(predic101[,6])
plot(reclama[,1]/reclama[,2],xlab="Categoría",ylab="Probabilidades",ylim=c(ymin,ymax))
points(predic102[,1],pch=4)
for (i  in 1:10){segments(i,predic102[i,4],i,predic102[i,6])}