{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Functions and Modeling Applications" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents:\n", "\n", "- [Functions and Modeling Applications](#Functions-and-Modeling-Applications) \n", " - [Creating Functions](#Creating-Functions) \n", " - [Loops](#Loops)\n", " - [Linear Algebra](#Linear-Algebra) \n", " - [Finite Markov Chains](#Finite-Markov-Chains) \n", "\n", "\n", "This lab covers: \n", "\n", "(1) User-defined functions;\n", "\n", "(2) Loops;\n", "\n", "(3) Linear algebra applications; \n", "\n", "(4) Modeling finite-state Markov chains.\n", "\n", "----" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating Functions " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this section we will cover the basics of creating custom functions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A function is essentially an object that takes a set of inputs, applies some sort of procedure to said inputs, and spits out a result.\n", "\n", "Functions can be handy for organizing code that is likely to be routinely re-used in the future with varying inputs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define a function named `add` that takes two variables, `x` and `y`, as inputs, and returns their sum." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "add (generic function with 1 method)" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function add(x, y)\n", " z = x + y\n", " return z\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To check, we call on `add()` with the inputs `x=2` and `y=3` -- obviously the output is supposed to be 5. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "add(2,3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Success!\n", "\n", "Recall that we stored the output into a variable `z` that was defined inside of the function.\n", "\n", "If we try to call on `z` in the global environment (outside of the function), we will get an error since `z` only lives within the function:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "ename": "LoadError", "evalue": "\u001b[91mUndefVarError: z not defined\u001b[39m", "output_type": "error", "traceback": [ "\u001b[91mUndefVarError: z not defined\u001b[39m", "", "Stacktrace:", " [1] top-level scope", " [2] include_string(::Function, ::Module, ::String, ::String) at .\\loading.jl:1091" ] } ], "source": [ "z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's define a function called `all_operations` that takes takes two variables, `x` and `y`, as inputs, and returns their sum, difference, product, and quotient. " ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "all_operations (generic function with 1 method)" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function all_operations(x, y)\n", " sum = x + y\n", " difference = x - y\n", " product = x * y\n", " quotient = x / y \n", " result = (sum, difference, product, quotient) \n", " return result \n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To check whether `all_operations` works, we call on it with the inputs `x=1` and `y=2`:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3, -1, 2, 0.5)" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "all_operations(1, 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the output of `all_operations()` is a tuple with four entires.\n", "\n", "Tuples are useful as function output objects because we can easily store their entries as separate variables:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "xy_sum = 3\n", "xy_difference = -1\n", "xy_product = 2\n", "xy_quotient = 0.5\n" ] } ], "source": [ "# Store output of `all_operations(1,2)` as separate variables\n", "xy_sum, xy_difference, xy_product, xy_quotient = all_operations(1,2)\n", "\n", "# Print all collected variables\n", "@show xy_sum\n", "@show xy_difference\n", "@show xy_product\n", "@show xy_quotient;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An alternative simpler way of defining the `all_operations` function by creating an equivalent `all_operations_v2`:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "all_operations_v2 (generic function with 1 method)" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function all_operations_v2(x,y)\n", " (x + y, x - y, x * y, x / y)\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What did we do differently?\n", "- We defined and stored all operation results in an unassigned tuple ;\n", "- We didn't use `return` at the end of the function to return the output.\n", "\n", "When we called on the function, Julia noticed the lack of a `return` command and chose the last thing it saw as the output -- in our case this was the unassigned tuple.\n", "\n", "Is this alternative way of defining `all_operations()` better? Not necessarily -- it depends on the context. For example, you might find that shorter code isn't *necessarily* easier to read.\n", "\n", "Let's just apply `all_operations()` and `all_operations_v2()` to the same inputs and check whether the outputs match using a custom-defined function `check()`: " ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"The two functions are the same!\"" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create var. `condition` that tests whether outputs are equivalent\n", "condition = all_operations(1,2) == all_operations_v2(1,2)\n", "\n", "# Create fun. `check` w/ input `condition`\n", "function check(condition)\n", " if condition == true \n", " result = \"The two functions are the same!\"\n", " end \n", " if condition != true\n", " result = \"The two functions are not the same!\"\n", " end \n", " return result\n", "end \n", "\n", "# Run `check` on `condition`\n", "check(condition)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `check()` function, as defined in the previous cell, is pretty clunky -- let's simplify it:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"The two functions are the same!\"" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Re-define function `check()`\n", "function check(condition)\n", " if condition == true\n", " return \"The two functions are the same!\"\n", " else \n", " return \"The two functions are not the same!\"\n", " end \n", "end\n", "\n", "# Call on `check()`\n", "check(condition)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or alternatively:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"The two functions are the same!\"" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Re-define function `check()\n", "function check(condition)\n", " if condition == true \n", " return \"The two functions are the same!\"\n", " end \n", " \"The two functions are not the same!\"\n", "end \n", "\n", "# Call on `check()`\n", "check(condition)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again -- in the case of simple functions such as the ones shown above, being super efficient is not necessary. \n", "\n", "But clunky code can make larger scripts hard to read, and potentially even run slow!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's talk math.\n", "\n", "Defining mathematical functions in Julia is easy.\n", "\n", "Let's define the polynomial mapping $f:\\mathbb{R} \\rightarrow \\mathbb{R}$ such that $f(x) = x^2 - 3x + 2$:" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "f (generic function with 1 method)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(x) = x^2 - 3x + 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we're interested in knowing the value of $f(\\pi)$:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.4448264403199786" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(pi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively:" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.4448264403199786" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f(π)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a lot of details on user-defined functions that we haven't been able to cover here, but I think the above should be enough to get you started. \n", "\n", "Visit the official Julia manual section on [functions](https://docs.julialang.org/en/v1/manual/functions/) to learn more." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loops" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's print every integer between 1 and 5 using a while loop:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n", "2\n", "3\n", "4\n", "5\n" ] } ], "source": [ "# Initial value\n", "i = 1\n", "\n", "# While loop:\n", "while i <= 5 # Run until i = 5\n", " println(i) # Print i\n", " i = i + 1 # Add 1 to i for the next iteration of the loop\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above can be accomplished more easily using a for loop:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n", "2\n", "3\n", "4\n", "5\n" ] } ], "source": [ "for i in 1:5 \n", " println(i)\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can pass any kind of sequence to a for loop.\n", "\n", "For example, we can print the set of odd numbers between 1 and 5 by defining a sequence called `sequence` and then iterating the values of said sequence:" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.0\n", "3.0\n", "5.0\n" ] } ], "source": [ "sequence = [1.0,3.0,5.0]\n", "\n", "for i in sequence \n", " println(i)\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we want to instead print the index associated with the entries of `sequence` instead of the entry values themselves?" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n", "2\n", "3\n" ] } ], "source": [ "for i in eachindex(sequence)\n", " println(i)\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose we want to square all values of `sequence` and store it as a separate array called `seq_out`.\n", "\n", "We can accomplish this using a for loop that goes through each entry of `sequence`, squares it, and stores it as the corresponding entry of `seq_out`:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 1.0\n", " 9.0\n", " 25.0" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Declare `seq_out` as a vector \n", "# w/ the same number of entries as `sequence`\n", "seq_out = zeros(length(sequence))\n", "\n", "# Run a for-loop that goes through\n", "# the indexes of `sequence`\n", "for i in eachindex(sequence) \n", " seq_out[i] = sequence[i]^2\n", "end \n", "\n", "# Print `seq_out`\n", "seq_out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, we can use a for loop in a **comprehension**:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 1.0\n", " 9.0\n", " 25.0" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "seq_out = [sequence[i]^2 for i in eachindex(sequence)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Even better -- we can **broadcast** (remember this from last lab?) `^2` across `sequence`." ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 1.0\n", " 9.0\n", " 25.0" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "seq_out = sequence.^2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear Algebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we load up the `LinearAlgebra` package." ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's assume we have vectors $a_1 = (1, 2, 3)'$ and $a_2 = (4, 5, 6)'$.\n", "\n", "We start by defining these two column vectors:" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "a_1 = [1; 2; 3]\n", "a_2 = [4, 5, 6];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Recall that whether we use `;`'s or `,`'s to separate entries in single-entry arrays, we still get column vectors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Obviously $a_1$ and $a_2$ do not span each other, but let's just check to be sure:" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " -0.0\n", " -0.0" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = [a_1 a_2]\n", "b = A \\ zeros(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the zero vector is the only solution for $x$ in $A \\, x = b$ where $A = [ a_1 a_2 ]$, then $a_1$ and $a_2$ must be linearly independent." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Find the dot product of $a_1$ and $a_2$:" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "32" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a_1' * a_2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively:" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "32" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dot(a_1, a_2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's find $a_1 a_2'$:" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Int64,2}:\n", " 4 5 6\n", " 8 10 12\n", " 12 15 18" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a_1 * a_2'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's add the two vectors:" ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 5\n", " 7\n", " 9" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a_1 + a_2 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Subtract $a_2$ from $a_1$:" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " -3\n", " -3\n", " -3" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a_1 - a_2 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's scale vector $a_1$ by 3:" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 3\n", " 6\n", " 9" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "3a_1 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equivalently:" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 3\n", " 6\n", " 9" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "3 * a_1 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equivalently:" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 3\n", " 6\n", " 9" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "3 .* a_1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The norm of vector $a_1$:" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.7416573867739413" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(a_1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since $a_1$ and $a_2$ cannot span $\\mathbb{R}^3$, we can find another orthogonal vector $a_3$.\n", "\n", "Is $a_3 = a_1 + a_2$ orthogonal? (Obviously not, by definition, but let's practice checking)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Array{Float64,1}:\n", " 0.9999999999999979\n", " 1.000000000000001" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a_3 = a_1 + a_2\n", "b = A \\ a_3 " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since there exists a non-trivial solution to $x$ in $[a_1 \\, a_2] \\, x = A \\, x = a_3$, then $a_3$ is not linearly independent. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can find a linearly independent $a_3$ by guessing some initial vector $b_3$, projecting it onto the columns of $A$ to obtain the projection $\\hat{b}_3$, and then extracting the orthogonal $a_3 = b_3 - \\hat{b}_3$:" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " -8.881784197001252e-16\n", " -1.3322676295501878e-15\n", " 2.6645352591003757e-15" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b_3 = [2, 3, 4]\n", "a_3 = b_3 - (A * inv(A'A) * A' * b_3)\n", "a_3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's check whether $a_3$ is actually linearly independent by redefining $A$ as $A = [a_1 \\, a_2 \\, a_3]$ and solving for $x$ in $A \\, x = 0$:" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0\n", " 0.0\n", " -0.0" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = [a_1 a_2 a_3]\n", "A \\ zeros(3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since the only solution for $x$ in $A \\, x = 0$ is the trivial solution, then $A$ must be full-rank." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's check the eigenvalues and eigenvectors of $A$:" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [], "source": [ "eigenv_A, eigenvec_A = eigen(A);" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " -0.4641016151377548\n", " 4.440892098500625e-15\n", " 6.464101615137752" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigenv_A" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3×3 Array{Float64,2}:\n", " -0.491831 2.96059e-16 -0.412884\n", " 0.180023 1.4803e-16 -0.56401\n", " 0.851877 1.0 -0.715136" ] }, "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigenvec_A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Finite Markov Chains" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose that for a Markov process with three states we are given an initial state distribution and a stochastic matrix. \n", "\n", "We are told to compute the state density in 10 periods (at $t=10$).\n", "\n", "The given initial distribution is $P_0 = (1/3, 1/3, 1/3)'$, while the stochastic matrix is \n", "\n", "$$M = \\begin{bmatrix} 0.95 & 0.05 & 0 \\\\ 0.15 & 0.75 & 0.1 \\\\ 0 & 0.5 & 0.5 \\end{bmatrix} \\, .$$ " ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.6415045988833332\n", " 0.2941181890520833\n", " 0.06437721206458333" ] }, "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define initial state distr.\n", "P0 = [1/3, 1/3, 1/3]\n", "\n", "# Define stochastic matrix\n", "M = [0.95 0.05 0 ; 0.15 0.75 0.1 ; 0.0 0.5 0.5]\n", "\n", "# Compute state distr. at t = 10\n", "P10 = (P0' * M^(10))'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More generally, what if we're interested in computing the state density for a variety of $n \\in \\mathbb{N}$ periods? \n", "\n", "This is when writing custom functions comes in handy!\n", "\n", "Let's create a function named `markov_chain` with inputs for $P_0$, $M$, and $n$:" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "markov_chain (generic function with 1 method)" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function markov_chain(P0, M, n)\n", " \n", " # Start by creating a copy of P0\n", " # to later feed into the loop\n", " P = copy(P0) \n", " \n", " # Run a loop that computes\n", " # P n-steps ahead\n", " for i in 1:n \n", " new_P = (P' * M)'\n", " P = new_P\n", " end \n", " \n", " # Return the final state distr.\n", " return P\n", "\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now try our new function out using the previously-defined `P0` and `M` with `n` set to 10:" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.641504598883333\n", " 0.29411818905208326\n", " 0.06437721206458331" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "markov_chain(P0, M, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Does this match `P10`?" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "true" ] }, "execution_count": 67, "metadata": {}, "output_type": "execute_result" } ], "source": [ "markov_chain(P0, M, 10) ≈ P10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It does!\n", "\n", "We can now use this function to do a whole bunch of useful things.\n", "\n", "For example, given the same $P_0$ and $M$, we can now create a list of $P_n$ for $n = 1,2,\\ldots,20$:" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [], "source": [ "probabilities = [markov_chain(P0, M, n) for n in 1:20];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We may visualize the probability of each state across time by gathering and plotting the corresponding series for all three states:" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [], "source": [ "state1 = zeros(20)\n", "state2 = zeros(20)\n", "state3 = zeros(20)\n", "\n", "for i in 1:20\n", " state1[i] = probabilities[i][1]\n", " state2[i] = probabilities[i][2]\n", " state3[i] = probabilities[i][3]\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's load the `Plots` package to make a couple of plots." ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [], "source": [ "using Plots" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's make a simple plot that contains all series:" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "time = 1:20\n", "plot(time, state1) # Plot state 1 density\n", "plot!(time, state2) # Plot state 2 density\n", "plot!(time, state3) # Plot state 3 density" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that with time, being in state 1 becomes more likely, while states 2 and 3 becomes less likely." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We may also create a separate plot for each series, but include them in a single composition:" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p1 = plot(time, state1, title = \"Normal Growth\") # Plot state 1 probability\n", "p2 = plot(time, state2, title = \"Recession\") # Plot state 2 probability\n", "p3 = plot(time, state3, title = \"Deep Recession\") # Plot state 3 probability\n", "plot(p1, p2, p3, layout = (3,1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's write a function that computes the probability of a given sequence of outcomes.\n", "\n", "(You should have seen something like this on Assignment 1.)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "outcome_prob (generic function with 1 method)" ] }, "execution_count": 81, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function outcome_prob(outcome, P0, M)\n", "\n", "# Make sure `outcome` contains integers\n", "# to allow for indexing \n", "outcome = floor.(Int64, outcome) \n", "\n", "# Store probability of initial state \n", "probability = P0[outcome[1]]\n", "\n", "# Compute probability of `outcome` sequence\n", "for i in 2:length(outcome)\n", " probability = probability * M[outcome[i-1], outcome[i]]\n", " end \n", " \n", " # Return `probability` -- prob. of sequence `outcome`\n", " return probability\n", "end " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's keep using our previously-defined initial distribution $P_0$ and stochastic matrix $M$.\n", "\n", "We can now feed the function the following sequence of states to check its functionality: $11$ -- the probability of this sequence of outcomes should obviously be (1/3)(19/20) according to our defined $P_0$ and $M$." ] }, { "cell_type": "code", "execution_count": 84, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "0.31666666666666665" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "outcome_prob(ones(2), P0, M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What if we want to find out the probability of having one of the following outcomes: (1) $11$, and (2) $22$ ?\n", "\n", "We can use `outcome_prob()` to compute the probability of each outcome, and then find their sum." ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.5666666666666667" ] }, "execution_count": 85, "metadata": {}, "output_type": "execute_result" } ], "source": [ "outcome_prob([1,1], P0, M) + outcome_prob([2,2], P0, M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.4", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.4" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }