Simulacion de Montecarlo

EJEMPLO

library(markovchain)

library(formattable)

#Definir la matriz P
estados = c(1:3)
P = matrix(c(0.5, 0.3, 0.2,
             0.3, 0.4, 0.3,
             0.2, 0.3, 0.5), nrow = 3, ncol = 3, byrow = T)

dimnames(P) = list(estados,estados)

#Estado actual
estado = 1

#Definir la distribución de probabilidad para el dia siguiente dado el estado actual
P[estado,]

#Obtener la calidad del aire para el dia siguiente
sample(estados, size = 1, prob = P[estado,])

#Simular 10 dias de la calidad del aire
dias = 10
set.seed(0)#Definiendo la semilla

#Definir el estado inicial
estado = 1

#Inicilizar una lista que contiene la calidad del aire en cada dia
lista_estados = c()

for (i in 1:dias) {
  #Guardar el estado actual en la lista de estados
  lista_estados = c(lista_estados,estado)
  
  #OBtener el estado futuro
  estado = sample(estados, size = 1, prob = P[estado,])
}
lista_estados


plot(lista_estados,type = 'l', xlab = 'Dias', ylab = 'Calidad del aire',
     main = 'Evolución de la calidad del aire')


# Simulación con 8 escenarios ---------------------------------------------


dias = 10
escenarios = 8
set.seed(0)

matriz_calidad = matrix(0,nrow = dias, ncol = escenarios)

for (j in 1:escenarios) {
  
  #Definir el estado actual
  estado = 1
  
  #Simular los 10 dias
  for (i in 1:dias) {
    
    #Guardar el estado en la matriz 
    matriz_calidad[i,j] = estado
    
    estado = sample(estados, size = 1, prob = P[estado,])
    
  }
}

matplot(matriz_calidad, type = 'l', xlab = 'Dias', ylab = 'Calidad del aire',
        main = 'Evolución de la calidad del aire')

EJERCICIO

#Tasa de la demanda
lambda<-15

#Modelación del manejo de inventarios de Wattenspharma con cadenas de Markov
#Crear los estados
estados<-c(0:100)

#*****Crear y llenar la matriz P de la política*****
matrizP<-matrix(0,nrow = length(estados), ncol = length(estados))
rownames(matrizP)<-estados
colnames(matrizP)<-estados

#Para la Politica -> si i<=70 solicita 30 cajas
for (i in estados) {
  for (j in estados) {
    if(i<=70 & j>0){
      matrizP[i+1,j+1]= dpois(30+i-j,lambda=lambda)
    }else if(i<=70 & j==0){
      matrizP[i+1,j+1]=ppois(30+i-1,lambda = lambda, lower.tail = F)
    }else if(i>70 & j>0){
      matrizP[i+1,j+1]=dpois(i-j,lambda=lambda)
    }else if(i>70 & j==0){
      matrizP[i+1,j+1]=ppois(i-1,lambda = lambda, lower.tail = F)
    }
  }
}

# Crear la cadena usando el paquete markovchain
cmtd<-new("markovchain", states=as.character(estados), transitionMatrix=matrizP)


# Literal B ---------------------------------------------------------------

# Vector para guardar el valor esperado del inventario en cada semana
vEsperadoInv <- c()    

# Recorrer las semanas
for(n in 1:30){
  estado_inicial <- "0"            # Definir el estado inicial
  P_t <- cmtd^n                    # Calcular la matriz de probabilidades en el transitorio
  probs <- P_t[estado_inicial,]    # Calcular probabilidades en el transitorio dado el estado inicial
  vEsperadoInv <- c(vEsperadoInv, probs%*%estados)  # Calcular valor esperado y agregar al vector
}

plot(vEsperadoInv, type = "l", ylab="Valor esperado del inventario", xlab="Semana", main="Inventario al final de la semana")

# Calcula el valor esperado del costo de ordenar de las siguientes 30 semanas
cOrdenar <- 5000
vEsperadoOrd <- 0

for(n in 1:30){
  estado_inicial <- "0"            # Definir el estado inicial
  P_t <- cmtd^n                    # Calcular la matriz de probabilidades en el transitorio
  probs <- P_t[estado_inicial,]    # Calcular probabilidades en el transitorio dado el estado inicial
  
  # Calcular el valor esperado de ordenar y sumarlo al costo hasta el momento
  vEsperadoOrd <- vEsperadoOrd + sum(probs[1:71])*cOrdenar
}
vEsperadoOrd



# Literal C. MonteCarlo ---------------------------------------------------
escenarios = 8
semanas = 30

inventario = matrix(0,nrow = semanas, ncol = escenarios)


#Simular la cadena de markov
set.seed(0)
for (j in 1:escenarios) {
  
  #Determinar el estado inicial
  estado = 0
  
  #SImular las 30 semanas
  for (i in 1:semanas) {
    
    #Guardar el inventario actual
    inventario[i,j] = estado
    
    estado = sample(estados, size = 1, prob = matrizP[as.character(estado),])
    
  }
}

matplot(inventario, type = 'l', xlab = 'Semanas', ylab = 'Inventario')


# Literal D ---------------------------------------------------------------


escenarios = 1000
semanas = 30

inventario = matrix(0,nrow = semanas, ncol = escenarios)


#Simular la cadena de markov
set.seed(0)
for (j in 1:escenarios) {
  
  #Determinar el estado inicial
  estado = 0
  
  #SImular las 30 semanas
  for (i in 1:semanas) {
    
    #Guardar el inventario actual
    inventario[i,j] = estado
    
    estado = sample(estados, size = 1, prob = matrizP[as.character(estado),])
    
  }
}

hist(inventario[30,], main = 'Distribución del inventario en la semana 30')

#Calculando el inventario promedio con MonteCarlo
inventario_promedio = mean(inventario[30,])

print(paste('Inventario con simulación de MonteCarlo: ', inventario_promedio))

#Valoresperado con analisis transitorio
print(paste('Valor esperado del inventario (literal b): ', round(vEsperadoInv[30],3)))



# Literal E ---------------------------------------------------------------


escenarios <- 1000
semanas <- 30

# Inicializar vector para guardar los costos de ordenar de cada escenario
vector_costos <- c()

# Simular la cadena de Markov
set.seed(0)
for(j in 1:escenarios){
  
  # Definir el estado inicial (inventario inicial)
  estado <- 0
  # Inicializar costo de ordenar
  costo <- 0
  
  # Simular las transiciones por 30 semanas
  for(i in 1:semanas){
    
    # Si el inventario es menor a 70, sumar costo de ordenar
    if(estado<=70){
      costo = costo + cOrdenar
    }
    estado <- sample(estados, size=1, prob=matrizP[as.character(estado),])
  }
  
  # Guardar el costo de ordenar
  vector_costos <- c(vector_costos, costo)
}
costo_promedio = mean(vector_costos)
costo_promedio

print(paste("Costo total de ordenar promedio utilizando simulación de Monte Carlo: ", currency(costo_promedio)))

print(paste("Valor esperado del costo total de ordenar: ", currency(round(vEsperadoOrd,3))))