# Cargar paquete
library(deSolve)
## Warning: package 'deSolve' was built under R version 4.4.3
# Parámetros
r <- 0.05 # tasa de interés
k0 <- 10 # capital inicial
T <- 2 # tiempo final
consumo_optimo <- 5.25 # consumo constante óptimo hallado
# Definir la ecuación diferencial: dk/dt = r*k - c
modelo <- function(t, estado, parametros) {
k <- estado[1]
c <- parametros["c"]
r <- parametros["r"]
dkdt <- r * k - c
return(list(c(dkdt)))
}
# Condiciones iniciales y parámetros
estado_inicial <- c(k = k0)
parametros <- c(c = consumo_optimo, r = r)
tiempo <- seq(0, T, by = 0.01)
# Resolver la ODE
resultado <- ode(y = estado_inicial, times = tiempo, func = modelo, parms = parametros)
# Convertir a data.frame para graficar
resultado_df <- as.data.frame(resultado)
print(resultado_df)
## time k
## 1 0.00 10.00000000
## 2 0.01 9.95248812
## 3 0.02 9.90495248
## 4 0.03 9.85739307
## 5 0.04 9.80980987
## 6 0.05 9.76220287
## 7 0.06 9.71457207
## 8 0.07 9.66691744
## 9 0.08 9.61923898
## 10 0.09 9.57153668
## 11 0.10 9.52381052
## 12 0.11 9.47606048
## 13 0.12 9.42828657
## 14 0.13 9.38048877
## 15 0.14 9.33266706
## 16 0.15 9.28482143
## 17 0.16 9.23695187
## 18 0.17 9.18905838
## 19 0.18 9.14114093
## 20 0.19 9.09319952
## 21 0.20 9.04523412
## 22 0.21 8.99724475
## 23 0.22 8.94923137
## 24 0.23 8.90119397
## 25 0.24 8.85313256
## 26 0.25 8.80504710
## 27 0.26 8.75693760
## 28 0.27 8.70880403
## 29 0.28 8.66064640
## 30 0.29 8.61246468
## 31 0.30 8.56425886
## 32 0.31 8.51602893
## 33 0.32 8.46777488
## 34 0.33 8.41949670
## 35 0.34 8.37119438
## 36 0.35 8.32286789
## 37 0.36 8.27451724
## 38 0.37 8.22614241
## 39 0.38 8.17774338
## 40 0.39 8.12932015
## 41 0.40 8.08087269
## 42 0.41 8.03240101
## 43 0.42 7.98390509
## 44 0.43 7.93538492
## 45 0.44 7.88684047
## 46 0.45 7.83827175
## 47 0.46 7.78967874
## 48 0.47 7.74106143
## 49 0.48 7.69241980
## 50 0.49 7.64375384
## 51 0.50 7.59506355
## 52 0.51 7.54634890
## 53 0.52 7.49760989
## 54 0.53 7.44884651
## 55 0.54 7.40005873
## 56 0.55 7.35124656
## 57 0.56 7.30240998
## 58 0.57 7.25354897
## 59 0.58 7.20466352
## 60 0.59 7.15575363
## 61 0.60 7.10681927
## 62 0.61 7.05786044
## 63 0.62 7.00887713
## 64 0.63 6.95986932
## 65 0.64 6.91083699
## 66 0.65 6.86178015
## 67 0.66 6.81269877
## 68 0.67 6.76359284
## 69 0.68 6.71446236
## 70 0.69 6.66530730
## 71 0.70 6.61612766
## 72 0.71 6.56692343
## 73 0.72 6.51769458
## 74 0.73 6.46844112
## 75 0.74 6.41916302
## 76 0.75 6.36986027
## 77 0.76 6.32053287
## 78 0.77 6.27118080
## 79 0.78 6.22180405
## 80 0.79 6.17240260
## 81 0.80 6.12297645
## 82 0.81 6.07352558
## 83 0.82 6.02404997
## 84 0.83 5.97454962
## 85 0.84 5.92502452
## 86 0.85 5.87547464
## 87 0.86 5.82589999
## 88 0.87 5.77630054
## 89 0.88 5.72667628
## 90 0.89 5.67702721
## 91 0.90 5.62735331
## 92 0.91 5.57765456
## 93 0.92 5.52793096
## 94 0.93 5.47818249
## 95 0.94 5.42840913
## 96 0.95 5.37861089
## 97 0.96 5.32878774
## 98 0.97 5.27893967
## 99 0.98 5.22906668
## 100 0.99 5.17916874
## 101 1.00 5.12924584
## 102 1.01 5.07929798
## 103 1.02 5.02932514
## 104 1.03 4.97932730
## 105 1.04 4.92930446
## 106 1.05 4.87925660
## 107 1.06 4.82918371
## 108 1.07 4.77908578
## 109 1.08 4.72896279
## 110 1.09 4.67881474
## 111 1.10 4.62864160
## 112 1.11 4.57844338
## 113 1.12 4.52822004
## 114 1.13 4.47797159
## 115 1.14 4.42769801
## 116 1.15 4.37739929
## 117 1.16 4.32707540
## 118 1.17 4.27672636
## 119 1.18 4.22635213
## 120 1.19 4.17595270
## 121 1.20 4.12552808
## 122 1.21 4.07507823
## 123 1.22 4.02460315
## 124 1.23 3.97410283
## 125 1.24 3.92357725
## 126 1.25 3.87302640
## 127 1.26 3.82245027
## 128 1.27 3.77184885
## 129 1.28 3.72122211
## 130 1.29 3.67057006
## 131 1.30 3.61989268
## 132 1.31 3.56918995
## 133 1.32 3.51846187
## 134 1.33 3.46770841
## 135 1.34 3.41692957
## 136 1.35 3.36612534
## 137 1.36 3.31529569
## 138 1.37 3.26444063
## 139 1.38 3.21356013
## 140 1.39 3.16265418
## 141 1.40 3.11172278
## 142 1.41 3.06076590
## 143 1.42 3.00978354
## 144 1.43 2.95877568
## 145 1.44 2.90774231
## 146 1.45 2.85668342
## 147 1.46 2.80559899
## 148 1.47 2.75448901
## 149 1.48 2.70335348
## 150 1.49 2.65219236
## 151 1.50 2.60100566
## 152 1.51 2.54979336
## 153 1.52 2.49855545
## 154 1.53 2.44729192
## 155 1.54 2.39600274
## 156 1.55 2.34468791
## 157 1.56 2.29334742
## 158 1.57 2.24198126
## 159 1.58 2.19058940
## 160 1.59 2.13917184
## 161 1.60 2.08772857
## 162 1.61 2.03625957
## 163 1.62 1.98476482
## 164 1.63 1.93324433
## 165 1.64 1.88169806
## 166 1.65 1.83012602
## 167 1.66 1.77852819
## 168 1.67 1.72690454
## 169 1.68 1.67525509
## 170 1.69 1.62357980
## 171 1.70 1.57187866
## 172 1.71 1.52015167
## 173 1.72 1.46839881
## 174 1.73 1.41662006
## 175 1.74 1.36481542
## 176 1.75 1.31298488
## 177 1.76 1.26112840
## 178 1.77 1.20924600
## 179 1.78 1.15733765
## 180 1.79 1.10540333
## 181 1.80 1.05344305
## 182 1.81 1.00145677
## 183 1.82 0.94944450
## 184 1.83 0.89740621
## 185 1.84 0.84534190
## 186 1.85 0.79325155
## 187 1.86 0.74113515
## 188 1.87 0.68899268
## 189 1.88 0.63682414
## 190 1.89 0.58462950
## 191 1.90 0.53240876
## 192 1.91 0.48016190
## 193 1.92 0.42788892
## 194 1.93 0.37558979
## 195 1.94 0.32326450
## 196 1.95 0.27091305
## 197 1.96 0.21853541
## 198 1.97 0.16613158
## 199 1.98 0.11370154
## 200 1.99 0.06124528
## 201 2.00 0.00876278
# Agregar columna de consumo constante
resultado_df$c <- consumo_optimo
# Graficar capital y consumo
plot(resultado_df$time, resultado_df$k, type = "l", col = "blue", lwd = 2,
xlab = "Tiempo (años)", ylab = "Capital / Consumo", main = "Optimización Dinámica")
lines(resultado_df$time, resultado_df$c, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Capital k(t)", "Consumo c(t)"),
col = c("blue", "red"), lty = c(1, 2), lwd = 2)
¿Qué hace el código?
Simula la evolución del capital \(k(t)\) en el tiempo.
Muestra que el consumo es constante.
Grafica cómo disminuye el capital para financiar ese consumo óptimo hasta llegar a cero al final del horizonte de tiempo.
-Objetivo: encontrar el nivel constante de consumo \(c\) que maximiza la utilidad total.
-Restricción: el capital evoluciona según una ecuación diferencial, dependiendo del consumo.
ecuacion <- expression(dk_dt == r * k - c)
print(ecuacion)
## expression(dk_dt == r * k - c)
library(deSolve)
# Parámetros
r <- 0.05
k0 <- 10
T <- 2
# Función para resolver la ODE dada c
resolver_capital <- function(c_val) {
modelo <- function(t, estado, parametros) {
k <- estado[1]
r <- parametros["r"]
c <- parametros["c"]
dkdt <- r * k - c
return(list(c(dkdt)))
}
estado_inicial <- c(k = k0)
parametros <- c(r = r, c = c_val)
tiempo <- seq(0, T, by = 0.01)
resultado <- ode(y = estado_inicial, times = tiempo, func = modelo, parms = parametros)
return(as.data.frame(resultado))
}
utilidad_total <- function(c_val) {
if (c_val <= 0) return(-Inf) # no se puede consumir negativo o cero
resultado_df <- resolver_capital(c_val)
# penalización si el capital cae a cero o negativo
if (any(resultado_df$k <= 0)) return(-1e6)
utilidad <- sum(log(c_val)) * 0.01 # utilidad con integración numérica (paso = 0.01)
return(utilidad)
}
opt_result <- optimize(utilidad_total, lower = 0.01, upper = r * k0, maximum = TRUE)
consumo_optimo <- opt_result$maximum
cat("Consumo óptimo:", consumo_optimo, "\n")
## Consumo óptimo: 0.4999559
resultado_final <- resolver_capital(consumo_optimo)
resultado_final$c <- consumo_optimo
plot(resultado_final$time, resultado_final$k, type = "l", col = "blue", lwd = 2,
xlab = "Tiempo", ylab = "Capital / Consumo", main = "Trayectoria Óptima")
lines(resultado_final$time, resultado_final$c, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Capital k(t)", "Consumo c(t)"),
col = c("blue", "red"), lty = c(1, 2), lwd = 2)
\(¿Que muestra la grafica?\)
Línea azul: evolución del capital k(t) con el consumo óptimo.
Línea roja punteada: el consumo constante óptimo c durante el tiempo.
El codigo muestra el valor óptimo de consumo constante,la trayectoria del capital bajo esa política y una visualización clara del problema de control óptimo.
Este código aplica una técnica de control óptimo numérico para encontrar el consumo constante óptimo que maximiza la utilidad total de un agente económico, sujeto a una restricción dinámica sobre el capital.