Euler_RK_method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal here is to show how we can implement Euler's method using Julia! \n",
    "Let's consider a disability insurance, this means that $S = \\{*, \\diamond, \\dagger\\}$, we want to find the survivalprobabilities\n",
    "\n",
    "$$\\begin{aligned}\n",
    "P(t,s) &= \n",
    "\\begin{bmatrix}\n",
    "    p_{**}(t,s)         & p_{*\\diamond}(t,s)        & p_{*\\dagger}(t,s) \\\\\n",
    "    p_{\\diamond*}(t,s)  & p_{\\diamond\\diamond}(t,s) & p_{\\diamond\\dagger}(t,s) \\\\\n",
    "    0                   &  0                        & 1\n",
    "\\end{bmatrix}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Λ (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#states: \n",
    "# 0: alive, 1: disabeld, 2: deceased\n",
    "\n",
    "# Transition rate matrix Λ\n",
    "function Λ(t)\n",
    "    #state0:\n",
    "    μ01(t) = 0.0004 + 10^(0.06*t-5.46)\n",
    "    μ02(t) = 0.0005 + 10^(0.038*t-4.12)\n",
    "    μ00(t) = -(μ01(t) + μ02(t))\n",
    "    #state1:\n",
    "    μ10(t) = 0.05 \n",
    "    μ12(t) = μ02(t)\n",
    "    μ11(t) = -(μ10(t)+μ12(t))\n",
    "    #state2:\n",
    "    # transition rates in the deceased state are zero\n",
    "    \n",
    "    L = [μ00(t) μ01(t) μ02(t)\n",
    "         μ10(t) μ11(t) μ12(t)\n",
    "         0       0     0    ]\n",
    "    \n",
    "    return L\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Euler (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function f(t,M)\n",
    "    return M*Λ(t)\n",
    "end\n",
    "\n",
    "function Euler(t0, n_states, h, tn)\n",
    "    \"\"\"\n",
    "    (t0): Int initial age \n",
    "    (n_states): Int, number of states\n",
    "    (h): Float, stepsize \n",
    "    (tn): Int, final age\n",
    "    \"\"\"\n",
    "    P0 = Matrix(1.0*I, n_states,n_states)\n",
    "\n",
    "    if t0 == tn\n",
    "        return P0\n",
    "    end\n",
    "\n",
    "    N = Int(round((tn-t0)/h))\n",
    "    D = Int(size(P0)[1])\n",
    "    # store N+1 (DxD)-matricies\n",
    "    P = zeros(Float64,D,D,N+1) \n",
    "    P[:,:,1] = P0 \n",
    "\n",
    "    \n",
    "    for n in 1:N\n",
    "        P[:,:, n+1] = P[:,:,n] + h*f(t0+n*h, P[:,:,n])\n",
    "    end\n",
    "\n",
    "    return P\n",
    "end "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0, 0.9998186582432559, 0.9996367084438109, 0.9994541437695489, 0.9992709573419496, 0.9990871422356756, 0.998902691478156, 0.9987175980491666, 0.9985318548804079, 0.998345454855079]"
     ]
    }
   ],
   "source": [
    "t0 = 30 \n",
    "n_states = 3\n",
    "h = 1/12     # stepsize \n",
    "tn = 120\n",
    "\n",
    "sol = Euler(t0, n_states, h, tn) # transition probabilities from 30 to 120 with stepsize h =1/12\n",
    "print(sol[1 , 1, 1:10])          # first 10 steps of survival\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RungeKutta method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "k4 (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function k1(t,M)\n",
    "    return f(t,M)\n",
    "end\n",
    "\n",
    "function k2(t,M)\n",
    "    return f(t+h/2, M +h*k1(t,M)/2)\n",
    "end\n",
    "\n",
    "function k3(t,M)\n",
    "    return f(t+h/2, M+ h*k2(t, M)/2)\n",
    "end \n",
    "\n",
    "function k4(t,M)\n",
    "    return f(t+h, M + h*k3(t,M))\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RK4 (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function RK4(t0, n_states, h, tn)\n",
    "    P0 = Matrix(1.0*I, n_states,n_states)\n",
    "\n",
    "    if t0 == P0 \n",
    "        return P0\n",
    "    end \n",
    "\n",
    "    N = Int(round((tn-t0)/h))\n",
    "    D = Int(size(P0)[1])\n",
    "    # store N+1 (DxD)-matricies\n",
    "    P = zeros(Float64,D,D,N+1) \n",
    "    P[:,:,1] = P0 \n",
    "\n",
    "    for n in 1:N\n",
    "        P[:,:,n+1] = P[:,:,n] + (h/6)*(k1(t0 + n*h, P[:,:,n]) + 2*k2(t0 +n*h, P[:,:,n]) +\n",
    "                                       2*k3(t0 +n*h, P[:,:, n]) + k4(t0 +n*h, P[:,:,n]))\n",
    "    end\n",
    "\n",
    "    return P\n",
    "end \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0, 0.999818354796958, 0.9996360981583859, 0.9994532232351914, 0.9992697231316997, 0.9990855909052406, 0.9989008195657305, 0.9987154020752523, 0.9985293313476322, 0.9983426002480129]"
     ]
    }
   ],
   "source": [
    "sol_RK = RK4(t0, n_states, h, tn)\n",
    "print(sol_RK[1,1,1:10]) "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.2",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}