---
title: "Paseo al Azar"
author: "Joaquín Ortega"
date: "Noviembre, 2018"
output: 
  pdf_document:
  fig_width: 6
fig_height: 4.5
lan: spanish
---

###  Guión para el curso de Probabilidad, Licenciatura en Matemáticas, UG. 
###  Tema 4: cadenas de Markov / Paseo al Azar


#############################################
#############################################
#############################################
#############################################
##  Cadenas de Markov
## Ejemplo del libro de Baclawski modificado.
### Paseo al azar con probabilidad de éxito 18/38 (ruleta americana), partiendo de 50 con barreras absorbentes en 0 y 100. Se grafica la distribución de probabilidad al paso n. La gráfica dinámica no se presenta en el pdf.

```{r, eval = FALSE}
p <- 18/38
q <- 1-p
n <- 100
c <- 50
m <- matrix(0,n+1,n+1)
for(i in 2:n){
  m[i,i-1] <- q
  m[i,i+1] <- p
}
m[1,1] <- 1
m[n+1,n+1] <- 1
x <- matrix(0,1,n+1)
x[1,c+1] <- 1
for(i in 1:600){
  x <- x%*% m
  plot((0:n)[x>0], x[x>0], xlim=c(0,n+1),lwd=8,
       ylim=c(0,0.5),xlab='i',ylab='P(X=i)',type='h',
       main='Distribución de probabilidad para un paseo al azar asimétrico')
  abline(v=50,col='red')
  text(25,0.4,paste('n=',i),cex=2)
  Sys.sleep(.1)
}
```

###  Trayectorias para el paseo al azar simétrico
```{r}
 set.seed(5678)
 nn <- 2000
 xi <- 2*rbinom(nn,1,0.5)-1
 sn <- cumsum(xi)
 plot(1:nn,sn,type='l',ylim=c(-2.5*sqrt(nn),2.5*sqrt(nn)),
      main='Paseo al azar simétrico',xlab='',ylab='')
 abline(h=0,col='red')
 lines(1:nn, sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
       type='l',col='red')
 lines(1:nn, -sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
       type='l',col='red')
 for(i in 1:50){
   lines(1:nn,cumsum(2*rbinom(nn,1,0.5)-1),type='l',
         col=rainbow(50)[i])
   Sys.sleep(.1)
 }
```

### Reducimos el número de trayectorias

```{r}
 set.seed(5678) 
 nn <- 2000 
 xi <- 2*rbinom(nn,1,0.5)-1 
 sn <- cumsum(xi) 
 plot(1:nn,sn,type='l',ylim=c(-2.5*sqrt(nn),2.5*sqrt(nn)), 
      main='Paseo al azar simétrico',xlab='',ylab='') 
 abline(h=0,col='red') 
 ceros <- numeric(5) 
 ceros[1] <- sum(sn== 0) 
 lines(1:nn, sqrt(2*(1:nn)*log(max(2,log(1:nn)))), 
       type='l',col='red') 
 lines(1:nn, -sqrt(2*(1:nn)*log(max(2,log(1:nn)))), 
       type='l',col='red') 
 for(i in 1:4){ 
   sn <- cumsum(2*rbinom(nn,1,0.5)-1) 
   lines(1:nn,sn,type='l',col=rainbow(5)[i]) 
   ceros[i+1] <- sum(sn==0) 
   Sys.sleep(.1)
 } 
 legend('topright',legend=ceros,col=c('black',
                      rainbow(5)[1:4]),lwd=2) 
``` 
  
## Ampliamos las trayectorias a 200,000 juegos 
```{r} 
 set.seed(5678) 
 nn <- 200000 
 xi <- 2*rbinom(nn,1,0.5)-1 
 sn <- cumsum(xi) 
 ceros <- numeric(5) 
 ceros[1] <- sum(sn== 0) 
 plot(1:nn,sn,type='l',ylim=c(-2.5*sqrt(nn),2.5*sqrt(nn)), 
      main='Paseo al azar simétrico',xlab='',ylab='') 
 abline(h=0,col='red') 
 lines(1:nn, sqrt(2*(1:nn)*log(max(2,log(1:nn)))), 
       type='l',col='red') 
 lines(1:nn, -sqrt(2*(1:nn)*log(max(2,log(1:nn)))), 
       type='l',col='red') 
 for(i in 1:4){ 
     sn <- cumsum(2*rbinom(nn,1,0.5)-1) 
     lines(1:nn,sn,type='l',col=rainbow(5)[i]) 
     ceros[i+1] <- sum(sn==0) 
     Sys.sleep(.05) 
 } 
 legend('topleft',legend=ceros,col=c('black',rainbow(5)[1:4]),
        lwd=2) 
``` 

### Ruina del jugador 
### p = 18/37 

###  Graficas de trayectorias 
```{r} 
set.seed(65432) 
nn <- 10000; pp <- 18/37 
sn <- cumsum(c(10,2*rbinom(nn,1,pp)-1)) 
nstop <- min(c(min((1:nn)[sn==0]),min((1:nn)[sn==20]),nn)) 
plot(0:(nstop-1), sn[1:nstop],ylim=c(0,20),pch=19, 
  xlim=c(0,200),type='l',xlab='',ylab='',main='Paseo al Azar')
 points(nstop-1,sn[nstop],pch=19) 
 abline(h=c(0,20),col='red') 
 abline(v=0,col='red') 
 for(i in 1:9){ 
  sn <- cumsum(c(10,2*rbinom(nn,1,pp)-1)) 
  nstop <- min(c(min((1:nn)[sn==0]),min((1:nn)[sn==20]),nn)) 
   lines(0:(nstop-1),sn[1:nstop],type='l',col=rainbow(35+i)[i]) 
   points(nstop-1,sn[nstop],pch=19,col=rainbow(50)[35+i]) 
 } 
``` 


### Aproximación probabilidad de ruina a=10, N=20 por simulación.

```{r} 
set.seed(1029) 
nncol = 1000; nn=2000; aa=10; NN=20; pp=18/37 
funcion0 <- function(x){min(match(0,x),nn,na.rm = TRUE)} 
funcionN <- function(x){min(match(NN,x),nn,na.rm = TRUE)} 
matt <- matrix(2*rbinom(nn*nncol,1,pp)-1,ncol=nncol) 
matt <- apply(matt,2,cumsum) 
matt <- rbind(rep(aa,nncol),matt +aa) 
vec0 <- apply(matt,2,funcion0)
vecN <- apply(matt,2,funcionN) 
sum(vec0<vecN) 
```

### Aprox. probabilidades de ruina para todos los a -->

```{r} 
set.seed(1029) 
nncol = 1000; nn=2000; NN=20; pp=18/37 
matt <- matrix(2*rbinom(nn*nncol,1,pp)-1,ncol=nncol) 
matt <- apply(matt,2,cumsum) 
probruina <- numeric(19) 
for(aa in 1:19){ 
mattaa <- rbind(rep(aa,nncol),matt +aa) 
vec0 <- apply(mattaa,2,funcion0) 
vecN <- apply(mattaa,2,funcionN) 
probruina[aa] <-  sum(vec0<vecN)/1000 
 } 

 plot(1:19,probruina,type='h',lwd=10,xlab='Capital inicial', 
      ylab='Probabilidad',ylim=c(0,1), 
      main='Probabilidad de Ruina en Ruleta Europea, N=20') 
  
 pruina.teo <- function(param,N){ 
   theta <- (1-param)/param 
   pruina <- numeric(N-1) 
   for(i in 1:(N-1)){ 
     pruina[i] <- (theta^N-theta^i)/(theta^N-1)} 
   return(pruina) 
 }     

 pruina20 <- pruina.teo(pp,20) 
 points(1:19,pruina20,pch=19,col='red') 
 abline(h=(0:10)/10,col='grey') 
``` 

###  Ley del arcoseno. Resultados Teóricos -->
### Ley discreta -->
```{r} 
pr0 <- function(k)  # Probabilidad de estar en 0 
{  choose(2*k,k)/2^{2*k} } 

asd <- function(n,k){ # Ley discreta del arcoseno 
 alpha=choose(2*k,k)*choose(2*(n-k),n-k)/2^{2*n}
 return(alpha) } 

 probs20 <- asd(10,0:9) 
 cumsum(probs20) 
``` 

 ###  Funcion de distribución continua 
```{r} 
arcseno <- function(x) {1/(pi*sqrt(x*(1-x)))} 
xval <- (1:100)/101 
plot(xval,arcseno(xval), ylim=c(0,2.5),type='l',lwd=2, 
      xlab='',ylab='',xaxp=c(0,1,10),main='Densidad de la  
      distribución arcoseno') 

 for(i in 0:10) 
   points(i/10,10*asd(10,i),pch=19,col='blue') 
 for(i in 0:20)
   points(i/20,20*asd(20,i),pch=19,col='red') 
 legend(0,0.3,c('n=10','n=20'),col=c('blue','red'), 
        pch=c(19,19)) 
``` 

### Función de distribución 
```{r} 
xval1 <- (0:1000)/1000 
plot(xval1,2*asin(sqrt(xval1))/pi,type='l',ylim=c(0,1),lwd=2, 
      xlab='',ylab='',main='Función de Distribución del Arcoseno') 
``` 

###  Ley del arcoseno. Simulación 
### Consideramos 20 juegos simétricos 
```{r}
 funcion.arc <- function(x){ 
   pos <- sum(x > 0) 
   ceros.arc <- which(x==0) 
   pos <-pos + sum(x[ceros.arc-1]>0) 
   return(pos) 
} 
set.seed(7584) 
 ll <- 20; nncol <- 100000 
 mat.arc <- matrix(2*rbinom(ll*nncol,1,0.5)-1,ncol=nncol)  
 mat.arc <- apply(mat.arc,2,cumsum) 
 vecpos <- numeric(nn) 
 vecpos <- apply(mat.arc,2,funcion.arc) 
 plot(table(vecpos)/nncol,lwd=6, xlab='',ylab='Probabilidad', 
      main='Ley del Arcoseno') 
 points((0:10)*2+.2,asd(10,0:10),col='red',pch=17,cex=1.5) 
``` 

###  La otra ley del arcoseno -->
###  Ultima visita al origen -->
```{r} 
set.seed(192837) 
nn = 20; nncol=100000 
matt <- matrix(2*rbinom(nn*nncol,1,0.5)-1,ncol=nncol) 
matt <- apply(matt,2,cumsum) 
matt <- rbind(rep(0,nncol),matt) 

funcion.max0 <- function(x){match(0,rev(x))} 
vec.arc <- apply(matt,2,funcion.max0) 
tab.arc <- (table(vec.arc)-1)/nncol 

plot(2*(0:(nn/2)), tab.arc,type='h',ylim=c(0,0.2), 
      lwd=4,ylab = '',xlab='', 
      main='Ley del arcoseno para la última visita') 
 points((0:10)*2+.2,asd(10,0:10),col='red',pch=17,cex=1.5) 
```