{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a072cd7",
   "metadata": {},
   "source": [
    "# Clase 2 — Análisis transitorio de CTMC\n",
    "# Regímenes de mercado financiero y laboral (decisión de portafolio)\n",
    "\n",
    "**Fuente.** Shi, J. (2026). *\"Optimal control of multiple Markov switching stochastic system with application to portfolio decision\"*.\n",
    "\n",
    "**Contexto del paper.** El paper modela a un inversionista que enfrenta **dos** mercados con cambio de régimen (\"regime switching\") gobernados por **dos cadenas de Markov distintas e independientes**: el mercado financiero $\\varepsilon(t)$ y el mercado laboral $\\zeta(t)$. En vez de tratar cada cadena por separado, el paper las combina en **una sola cadena conjunta** $\\xi(t)$ con 4 estados, y deriva explícitamente su matriz generadora combinada (su Ecuación 3.2). Hoy reproducimos exactamente esa construcción con datos ilustrativos, y la usamos para calcular **probabilidades transientes**: ¿qué tan probable es, en el año $t$, que la economía esté en cada combinación de régimen financiero/laboral?\n",
    "\n",
    "Curso: Modelado de Sistemas bajo Incertidumbre · Universidad de los Andes · 2026-20\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dcf79cc",
   "metadata": {},
   "source": [
    "## 1. Descripción del sistema\n",
    "\n",
    "Un inversionista observa dos mercados que cambian de régimen de forma continua e independiente entre sí:\n",
    "\n",
    "- **Mercado financiero:** alterna entre un régimen \"Alcista\" (Bull) y uno \"Bajista\" (Bear).\n",
    "- **Mercado laboral:** alterna entre \"Expansión\" (creación de empleo) y \"Recesión\" (destrucción de empleo).\n",
    "\n",
    "**Dos variables de estado:**\n",
    "- $\\varepsilon(t)\\in\\{0,1\\}$: régimen financiero (0=Alcista, 1=Bajista).\n",
    "- $\\zeta(t)\\in\\{0,1\\}$: régimen laboral (0=Expansión, 1=Recesión).\n",
    "\n",
    "## 2. Espacio de estados\n",
    "\n",
    "Siguiendo la notación del paper, el estado conjunto se numera como $\\xi(t)\\in\\{1,2,3,4\\}$:\n",
    "\n",
    "$$\\xi=1 \\text{ si } (\\varepsilon,\\zeta)=(0,0) \\qquad \\xi=2 \\text{ si } (\\varepsilon,\\zeta)=(1,0) \\qquad \\xi=3 \\text{ si } (\\varepsilon,\\zeta)=(0,1) \\qquad \\xi=4 \\text{ si } (\\varepsilon,\\zeta)=(1,1)$$\n",
    "\n",
    "## 3. Datos / parámetros\n",
    "\n",
    "Duraciones medias típicas de cada régimen (valores ilustrativos, del orden de magnitud de ciclos económicos reales):\n",
    "\n",
    "**Mercado financiero ($\\varepsilon$):**\n",
    "- Alcista: duración media 4 años \n",
    "- Bajista: duración media 1.5 años \n",
    "\n",
    "**Mercado laboral ($\\zeta$):**\n",
    "- Expansión: duración media 6 años \n",
    "- Recesión: duración media 2 años "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79ea24f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# importar librerías\n",
    "import numpy as np\n",
    "from jmarkov.ctmc import ctmc\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5737b53a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 1: define las tasas de cada mercado\n",
    "lda0_eps, lda1_eps = None, None\n",
    "lda0_zeta, lda1_zeta = None, None\n",
    "print(f'financiero: lda0={lda0_eps:.4f}, lda1={lda1_eps:.4f}')\n",
    "print(f'laboral:    lda0={lda0_zeta:.4f}, lda1={lda1_zeta:.4f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b800fa67",
   "metadata": {},
   "source": [
    "## 4. Vector de estados\n",
    "\n",
    "Ordenando $(\\varepsilon,\\zeta)$ igual que el paper (варía $\\varepsilon$ primero):\n",
    "\n",
    "$$\\text{estados}=[(0,0),(1,0),(0,1),(1,1)] \\;\\equiv\\; [\\xi{=}1,\\ \\xi{=}2,\\ \\xi{=}3,\\ \\xi{=}4]$$\n",
    "\n",
    "## 5. Vector de condiciones iniciales\n",
    "\n",
    "La economía arranca en Alcista + Expansión: $\\alpha=(1,0,0,0)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f5ac019",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 2: vector de estados y alpha\n",
    "estados = [(e, z) for z in range(2) for e in range(2)]\n",
    "idx = {s: i for i, s in enumerate(estados)}\n",
    "n = len(estados)\n",
    "\n",
    "alpha = np.zeros(n)\n",
    "alpha[idx[(0, 0)]] = None\n",
    "print(f'estados={estados}')\n",
    "print(f'alpha={alpha}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6318b28",
   "metadata": {},
   "source": [
    "## 6. Gráfico de tasas\n",
    "\n",
    "```\n",
    "              lda0_eps\n",
    "   (0,0) -----------------> (1,0)      [ Alcista+Expansion <-> Bajista+Expansion ]\n",
    "     |  <----------------- lda1_eps          |\n",
    "lda0_zeta | lda1_zeta               lda0_zeta | lda1_zeta\n",
    "     v                                          v\n",
    "   (0,1) -----------------> (1,1)      [ Alcista+Recesion  <-> Bajista+Recesion  ]\n",
    "              lda0_eps\n",
    "              <----------------- lda1_eps\n",
    "```\n",
    "\n",
    "Los movimientos horizontales cambian el régimen financiero ($\\varepsilon$); los verticales cambian el régimen laboral ($\\zeta$). Nunca cambian los dos al mismo tiempo, porque las cadenas son independientes y continuas.\n",
    "\n",
    "## 7. Tasas de transición\n",
    "\n",
    "$$Q^\\varepsilon=\\begin{pmatrix} \\end{pmatrix} \\qquad Q^\\zeta=\\begin{pmatrix}\\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b3f6cb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 3: construye Q_eps (2x2) y Q_zeta (2x2)\n",
    "Q_eps = np.array([\n",
    "\n",
    "])\n",
    "Q_zeta = np.array([\n",
    "\n",
    "])\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(Q_eps.sum(axis=1), 0) and np.allclose(Q_zeta.sum(axis=1), 0)\n",
    "print('OK: Q_eps y Q_zeta son matrices generadoras válidas.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93a27ba5",
   "metadata": {},
   "source": [
    "## 8. Matriz generadora conjunta\n",
    "\n",
    "Como el paper asume que $\\varepsilon(t)$ y $\\zeta(t)$ son independientes, el generador conjunto es la suma de Kronecker (equivalente, en este caso independiente, a la Ecuación 3.2 del paper):\n",
    "\n",
    "$$Q=Q^\\zeta\\otimes I_\\varepsilon + I_\\zeta\\otimes Q^\\varepsilon$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "613c7968",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 4: construye Q con la suma de Kronecker (orden: eps varía más rápido)\n",
    "I_e, I_z = np.eye(2), np.eye(2)\n",
    "Q = None  # usa np.kron(Q_zeta, I_e) + np.kron(I_z, Q_eps)\n",
    "print(f'Q=\\n{np.round(Q,4)}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(Q.sum(axis=1), 0)\n",
    "print('OK: Q es una matriz generadora válida (4x4).')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1de56381",
   "metadata": {},
   "source": [
    "## 9. Tasa total de salida de cada estado\n",
    "\n",
    "$q_i=-Q_{ii}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "713bc8f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 5: tasa total de salida de cada estado\n",
    "q = None  # usa -np.diag(Q)\n",
    "for s, qi in zip(estados, q):\n",
    "    print(f'estado {s}: q = {qi:.4f} por año')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdbb9ecc",
   "metadata": {},
   "source": [
    "## 10. Tiempos de permanencia\n",
    "\n",
    "$E[T_i]=1/q_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "023d077b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 6: tiempo medio de permanencia en cada estado conjunto\n",
    "E_T = None  # usa 1/q\n",
    "for s, t in zip(estados, E_T):\n",
    "    print(f'estado {s}: E[T] = {t:.4f} años')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "919d2904",
   "metadata": {},
   "source": [
    "## 11. Probabilidades de la cadena embebida\n",
    "\n",
    "$P_{ij}=q_{ij}/q_i$ para $i\\neq j$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90fe68b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 7: matriz de la cadena embebida\n",
    "P = np.zeros_like(Q)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        if i != j:\n",
    "            P[i, j] = None  # Q[i,j] / q[i]\n",
    "print(f'P (cadena embebida)=\\n{np.round(P,4)}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.allclose(P.sum(axis=1), 1)\n",
    "print('OK: cada fila de P suma 1.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0315a4c5",
   "metadata": {},
   "source": [
    "## 12. Probabilidades transientes $\\pi(t)$\n",
    "\n",
    "Con `jmarkov`, evaluamos $\\pi(t)=\\alpha\\,e^{Qt}$ para varios horizontes de inversión $t$ (en años)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b991ef69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 8: construye el modelo y evalúa pi(t)\n",
    "mc = ... # Construya la CTMC con la matriz generadora Q\n",
    "tiempos = [1, 3, 5, 10]\n",
    "for t in tiempos:\n",
    "    pi_t = None  # usa mc.transient_probabilities(t, alpha)\n",
    "    print(f't={t} años:')\n",
    "    for s, p in zip(estados, pi_t):\n",
    "        print(f'   pi{s}({t}) = {p:.5f}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51137957",
   "metadata": {},
   "source": [
    "## 13. Probabilidades marginales y peor escenario conjunto\n",
    "\n",
    "La probabilidad marginal de cada mercado se obtiene sumando $\\pi(t)$ sobre la otra variable:\n",
    "\n",
    "$$P(\\varepsilon(t)=e)=\\sum_{z}\\pi_{(e,z)}(t) \\qquad P(\\zeta(t)=z)=\\sum_{e}\\pi_{(e,z)}(t)$$\n",
    "\n",
    "El \"peor escenario\" para el inversionista es el estado conjunto Bajista + Recesión, $(\\varepsilon,\\zeta)=(1,1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf19c7ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 9: probabilidades marginales y probabilidad del peor escenario en t=10\n",
    "t = 10\n",
    "pi_t = ... # calcule las pribabilidades transientes\n",
    "\n",
    "P_eps = {e: None for e in range(2)}   # suma pi_t[idx[(e,z)]] sobre z\n",
    "P_zeta = {z: None for z in range(2)}  # suma pi_t[idx[(e,z)]] sobre e\n",
    "print(f'P(financiero(t)=e) en t={t}: {P_eps}')\n",
    "print(f'P(laboral(t)=z) en t={t}: {P_zeta}')\n",
    "\n",
    "P_peor = None  # pi_t en el estado (1,1)\n",
    "print(f'P(Bajista y Recesion) en t={t} = {P_peor:.5f}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.isclose(sum(P_eps.values()), 1) and np.isclose(sum(P_zeta.values()), 1)\n",
    "print('OK: las probabilidades marginales suman 1.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e2d5cc7",
   "metadata": {},
   "source": [
    "## 14. Probabilidad del peor escenario en el tiempo (gráfico)\n",
    "\n",
    "Evalúa $P(\\varepsilon(t)=1,\\zeta(t)=1)$ en una malla fina de tiempos y grafica cómo converge a su valor estacionario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "374e5073",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 10: evalúa P(peor escenario)(t) en una malla fina y grafica\n",
    "t_fino = np.linspace(0, 20, 150)\n",
    "P_peor_fino = None  # arreglo con P(eps=1,zeta=1)(t) para cada t en t_fino\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4.5))\n",
    "ax.plot(t_fino, P_peor_fino, color='#4C72B0', lw=2)\n",
    "ax.axhline(P_peor_fino[-1], color='#C44E52', ls='--', lw=1, label=f'Valor estacionario = {P_peor_fino[-1]:.4f}')\n",
    "ax.set_xlabel('t (años)')\n",
    "ax.set_ylabel('P(Bajista y Recesion)(t)')\n",
    "ax.set_title('Probabilidad del peor escenario conjunto vs. tiempo')\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b96adcc",
   "metadata": {},
   "source": [
    "## 15. Distribución estacionaria y medidas de desempeño"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15523381",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 11: distribución estacionaria conjunta y verificación por independencia\n",
    "ss = None  # usa mc.steady_state()\n",
    "print(f'pi (orden {estados}) = {np.round(ss,5)}')\n",
    "\n",
    "P_peor_ss = None\n",
    "print(f'P(Bajista y Recesion) estacionaria = {P_peor_ss:.5f}')\n",
    "\n",
    "# verificación: al ser independientes, P(eps=1,zeta=1) = P(eps=1)*P(zeta=1) en estado estacionario\n",
    "P_bear_ss = None  # usa ctmc(Q_eps).steady_state()\n",
    "P_rec_ss = None   # usa ctmc(Q_zeta).steady_state()\n",
    "print(f'P(Bajista)={P_bear_ss:.5f}  P(Recesion)={P_rec_ss:.5f}  producto={P_bear_ss*P_rec_ss:.5f}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert np.isclose(P_peor_ss, P_bear_ss * P_rec_ss, atol=1e-6)\n",
    "print('OK: la probabilidad conjunta estacionaria se factoriza (independencia).')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d74962dc",
   "metadata": {},
   "source": [
    "## 16. Costeo: costo esperado a partir de probabilidades\n",
    "\n",
    "Supongamos que cada combinación de régimen tiene asociado un **costo** (o pérdida esperada) para el inversionista, por ejemplo en miles de USD por año. Definimos una tasa de costo $c_{(e,z)}$ para cada estado conjunto:\n",
    "\n",
    "| Estado $(\\varepsilon,\\zeta)$ | Costo $c$ |\n",
    "|---|---|\n",
    "| (0,0) Alcista+Expansión | 0 |\n",
    "| (1,0) Bajista+Expansión | 2 |\n",
    "| (0,1) Alcista+Recesión | 1 |\n",
    "| (1,1) Bajista+Recesión | 5 |\n",
    "\n",
    "Nótese que $c_{(1,1)}=5 \\neq c_{(1,0)}+c_{(0,1)}=2+1=3$: el costo **no es aditivo**, hay un costo adicional por la interacción de ambas crisis a la vez. Esto es exactamente lo que distingue calcular el costo con la **distribución conjunta** de intentar aproximarlo solo con las **marginales**.\n",
    "\n",
    "**Costo esperado con la distribución conjunta (correcto):**\n",
    "\n",
    "$$E[C(t)]=\\sum_{(e,z)} c_{(e,z)}\\,\\pi_{(e,z)}(t)$$\n",
    "\n",
    "**Aproximación aditiva usando solo las marginales (válida solo si el costo fuera aditivo):**\n",
    "\n",
    "$$\\hat{E}[C(t)]=\\sum_e c^{(1)}_e\\,P(\\varepsilon(t)=e) + \\sum_z c^{(2)}_z\\,P(\\zeta(t)=z), \\qquad c^{(1)}=(0,2),\\ c^{(2)}=(0,1)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6c2272a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO 12: costo esperado conjunto vs. aproximación aditiva con marginales, para varios t\n",
    "c_joint = {(0, 0): 0, (1, 0): 2, (0, 1): 1, (1, 1): 5}\n",
    "c1 = {0: 0, 1: 2}   # costo marginal atribuible al régimen financiero\n",
    "c2 = {0: 0, 1: 1}   # costo marginal atribuible al régimen laboral\n",
    "\n",
    "for t in [1, 3, 5, 10]:\n",
    "    pi_t = mc.transient_probabilities(t, alpha)\n",
    "    E_C_joint = None  # suma de c_joint[s] * pi_t[idx[s]] sobre todos los estados s\n",
    "\n",
    "    P_e = {e: sum(pi_t[idx[(e, z)]] for z in range(2)) for e in range(2)}\n",
    "    P_z = {z: sum(pi_t[idx[(e, z)]] for e in range(2)) for z in range(2)}\n",
    "    E_C_aditivo = None  # suma c1[e]*P_e[e] + suma c2[z]*P_z[z]\n",
    "\n",
    "    print(f't={t}: E[C(t)] conjunto={E_C_joint:.4f}   aproximacion aditiva={E_C_aditivo:.4f}   diferencia={E_C_joint-E_C_aditivo:.4f}')\n",
    "\n",
    "# costo esperado de largo plazo (estado estacionario)\n",
    "E_C_ss = None  # suma de c_joint[s] * ss[idx[s]] sobre todos los estados s\n",
    "print(f'\\nCosto esperado de largo plazo (estacionario) = {E_C_ss:.4f}')\n",
    "\n",
    "# --- prueba rápida (no la modifiques) ---\n",
    "assert E_C_joint > E_C_aditivo, 'el costo conjunto debe superar la aproximacion aditiva (hay interaccion positiva)'\n",
    "print('OK: el costo conjunto supera la aproximacion aditiva, por el termino de interaccion en (1,1).')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff9d42e",
   "metadata": {},
   "source": [
    "## Conclusión\n",
    "\n",
    "_¿Por qué la probabilidad conjunta estacionaria del peor escenario es menor que la probabilidad de \"solo Bajista\" o \"solo Recesión\" por separado? Si el inversionista solo pudiera observar el mercado financiero (no el laboral), ¿qué información perdería?_"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}