getting_started_julia/julia_by_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='julia-by-example'></a>\n",
    "<div id=\"qe-notebook-header\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introductory Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Introductory Examples](#Introductory-Examples)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Example: Plotting a White Noise Process](#Example:-Plotting-a-White-Noise-Process)  \n",
    "  - [Example: Variations on Fixed Points](#Example:-Variations-on-Fixed-Points)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "We’re now ready to start learning the Julia language itself"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Level\n",
    "\n",
    "Our approach is aimed at those who already have at least some knowledge of programming — perhaps experience with Python, MATLAB, Fortran, C or similar\n",
    "\n",
    "In particular, we assume you have some familiarity with fundamental programming concepts such as\n",
    "\n",
    "- variables  \n",
    "- arrays or vectors  \n",
    "- loops  \n",
    "- conditionals (if/else)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approach\n",
    "\n",
    "In this lecture we will write and then pick apart small Julia programs\n",
    "\n",
    "At this stage the objective is to introduce you to basic syntax and data structures\n",
    "\n",
    "Deeper concepts—how things work—will be covered in later lectures\n",
    "\n",
    "Since we are looking for simplicity the examples are a little contrived\n",
    "\n",
    "In this lecture, we will often start with a direct MATLAB/FORTRAN approach which often is **poor coding style** in Julia, but then move towards more **elegant code** which is tightly connected to the mathematics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set Up\n",
    "\n",
    "We assume that you’ve worked your way through [our getting started lecture](getting_started.ipynb) already\n",
    "\n",
    "In particular, the easiest way to install and precompile all the Julia packages used in QuantEcon\n",
    "notes is to type `] add InstantiateFromURL` and then work in a Jupyter notebook, as described [here](getting_started.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other References\n",
    "\n",
    "The definitive reference is [Julia’s own documentation](https://docs.julialang.org/en/v1/)\n",
    "\n",
    "The manual is thoughtfully written but is also quite dense (and somewhat evangelical)\n",
    "\n",
    "The presentation in this and our remaining lectures is more of a tutorial style based around examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Plotting a White Noise Process\n",
    "\n",
    "To begin, let’s suppose that we want to simulate and plot the white noise\n",
    "process $ \\epsilon_0, \\epsilon_1, \\ldots, \\epsilon_T $, where each draw $ \\epsilon_t $ is independent standard normal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction to Packages\n",
    "\n",
    "The first step is to activate a project environment, which is encapsulated by `Project.toml` and `Manifest.toml` files\n",
    "\n",
    "There are three ways to install packages and versions (where the first two methods are discouraged, since they may lead to package versions out-of-sync with the notes)\n",
    "\n",
    "1. `add` the packages directly into your global installation (e.g. `Pkg.add(\"MyPackage\")` or `] add MyPackage`)  \n",
    "1. download an `Project.toml` and `Manifest.toml` file in the same directory as the notebook (i.e. from the `@__DIR__` argument), and then call `using Pkg; Pkg.activate(@__DIR__);`  \n",
    "1. use the `InstantiateFromURL` package  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "using InstantiateFromURL\n",
    "github_project(\"QuantEcon/quantecon-notebooks-julia\", version = \"0.2.0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have never run this code on a particular computer, it is likely to take a long time as it downloads, installs, and compiles all dependent packages\n",
    "\n",
    "This code will download and install project files from GitHub, [QuantEcon/QuantEconLecturePackages](https://github.com/QuantEcon/QuantEconLecturePackages/)\n",
    "\n",
    "We will discuss it more in [Tools and Editors](more_julia/tools_editors.ipynb#tools-editors), but these files provide a listing of packages and versions used by the code\n",
    "\n",
    "This ensures that an environment for running code is **reproducible**, so that anyone can replicate the precise set of package and versions used in construction\n",
    "\n",
    "The careful selection of package versions is crucial for reproducibility, as otherwise your code can be broken by changes to packages out of your control\n",
    "\n",
    "After the installation and activation, `using` provides a way to say that a particular code or notebook will use the package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using LinearAlgebra, Statistics, Compat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='import'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using Functions from a Package\n",
    "\n",
    "Some functions are built into the base Julia, such as `randn`, which returns a single draw from a normal distibution with mean 0 and variance 1 if given no parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "randn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other functions require importing all of the names from an external library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using Plots\n",
    "gr(fmt=:png); # setting for easier display in jupyter notebooks\n",
    "\n",
    "n = 100\n",
    "ϵ = randn(n)\n",
    "plot(1:n, ϵ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s break this down and see how it works\n",
    "\n",
    "The effect of the statement `using Plots` is to make all the names exported by the `Plots` module available\n",
    "\n",
    "Because we used `Pkg.activate` previously, it will use whatever version of `Plots.jl` that was specified in the `Project.toml` and `Manifest.toml` files\n",
    "\n",
    "The other packages `LinearAlgebra` and `Statistics` are base Julia libraries, but require an explicit using\n",
    "\n",
    "The arguments to `plot` are the numbers `1,2, ..., n` for the x-axis, a vector `ϵ` for the y-axis, and (optional) settings\n",
    "\n",
    "The function `randn(n)` returns a column vector `n` random draws from a normal distribution with mean 0 and variance 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arrays\n",
    "\n",
    "As a language intended for mathematical and scientific computing, Julia has\n",
    "strong support for using unicode characters\n",
    "\n",
    "In the above case, the `ϵ` and many other symbols can be typed in most Julia editor by providing the LaTeX and `<TAB>`, i.e. `\\epsilon<TAB>`\n",
    "\n",
    "The return type is one of the most fundamental Julia data types: an array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "typeof(ϵ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϵ[1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The information from `typeof()` tells us that `ϵ` is an array of 64 bit floating point values, of dimension 1\n",
    "\n",
    "In Julia, one-dimensional arrays are interpreted as column vectors for purposes of linear algebra\n",
    "\n",
    "The `ϵ[1:5]` returns an array of the first 5 elements of `ϵ`\n",
    "\n",
    "Notice from the above that\n",
    "\n",
    "- array indices start at 1 (like MATLAB and Fortran, but unlike Python and C)  \n",
    "- array elements are referenced using square brackets (unlike MATLAB and Fortran)  \n",
    "\n",
    "\n",
    "To get **help and examples** in Jupyter or other julia editor, use the `?` before a function name or syntax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```julia\n",
    "?typeof\n",
    "\n",
    "search: typeof typejoin TypeError\n",
    "\n",
    "Get the concrete type of x.\n",
    "\n",
    "Examples\n",
    "\n",
    "julia> a = 1//2;\n",
    "\n",
    "julia> typeof(a)\n",
    "Rational{Int64}\n",
    "\n",
    "julia> M = [1 2; 3.5 4];\n",
    "\n",
    "julia> typeof(M)\n",
    "Array{Float64,2}\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For Loops\n",
    "\n",
    "Although there’s no need in terms of what we wanted to achieve with our\n",
    "program, for the sake of learning syntax let’s rewrite our program to use a\n",
    "`for` loop for generating the data\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">In Julia v0.7 and up, the rules for variables accessed in `for` and `while` loops can be sensitive to how they are used (and variables can sometimes require a `global` as part of the declaration).  We strongly advise you to avoid top level (i.e. in the REPL or outside of functions) `for` and `while` loops outside of Jupyter notebooks.  This issue does not apply when used within functions\n",
    "\n",
    "Starting with the most direct version, and pretending we are in a world where `randn` can only return a single value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# poor style\n",
    "n = 100\n",
    "ϵ = zeros(n)\n",
    "for i in 1:n\n",
    "    ϵ[i] = randn()\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we first declared `ϵ` to be a vector of `n` numbers, initialized by the floating point `0.0`\n",
    "\n",
    "The `for` loop then populates this array by successive calls to `randn()`\n",
    "\n",
    "Like all code blocks in Julia, the end of the `for` loop code block (which is just one line here) is indicated by the keyword `end`\n",
    "\n",
    "The word `in` from the `for` loop can be replaced by either `∈` or `=`\n",
    "\n",
    "The index variable is looped over for all integers from `1:n` – but this does not actually create a vector of those indices\n",
    "\n",
    "Instead, it creates an **iterator** that is looped over – in this case the **range** of integers from `1` to `n`\n",
    "\n",
    "While this example successfully fills in `ϵ` with the correct values, it is very indirect as the connection between the index `i` and the `ϵ` vector is unclear\n",
    "\n",
    "To fix this, use `eachindex`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# better style\n",
    "n = 100\n",
    "ϵ = zeros(n)\n",
    "for i in eachindex(ϵ)\n",
    "    ϵ[i] = randn()\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, `eachindex(ϵ)` returns an iterator of indices which can be used to access `ϵ`\n",
    "\n",
    "While iterators are memory efficient because the elements are generated on the fly rather than stored in memory, the main benefit is (1) it can lead to code which is clearer and less prone to typos; and (2) it allows the compiler flexibility to creatively generate fast code\n",
    "\n",
    "In Julia you can also loop directly over arrays themselves, like so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϵ_sum = 0.0 # careful to use 0.0 here, instead of 0\n",
    "m = 5\n",
    "for ϵ_val in ϵ[1:m]\n",
    "    ϵ_sum = ϵ_sum + ϵ_val\n",
    "end\n",
    "ϵ_mean = ϵ_sum / m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where `ϵ[1:m]` returns the elements of the vector at indices `1` to `m`\n",
    "\n",
    "Of course, in Julia there are built in functions to perform this calculation which we can compare against"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϵ_mean ≈ mean(ϵ[1:m])\n",
    "ϵ_mean ≈ sum(ϵ[1:m]) / m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In these examples, note the use of `≈` to test equality, rather than `==`, which is appropriate for integers and other types\n",
    "\n",
    "Approximately equal, typed with `\\approx<TAB>`, is the appropriate way to compare any floating point numbers due to the standard issues of [floating point math](https://floating-point-gui.de/)\n",
    "\n",
    "\n",
    "<a id='user-defined-functions'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### User-Defined Functions\n",
    "\n",
    "For the sake of the exercise, let’s go back to the `for` loop but restructure our program so that generation of random variables takes place within a user-defined function\n",
    "\n",
    "To make things more interesting, instead of directly plotting the draws from the distribution, let’s plot the squares of these draws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# poor style\n",
    "function generatedata(n)\n",
    "    ϵ = zeros(n)\n",
    "    for i in eachindex(ϵ)\n",
    "        ϵ[i] = (randn())^2 # squaring the result\n",
    "    end\n",
    "    return ϵ\n",
    "end\n",
    "\n",
    "data = generatedata(10)\n",
    "plot(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here\n",
    "\n",
    "- `function` is a Julia keyword that indicates the start of a function definition  \n",
    "- `generatedata` is an arbitrary name for the function  \n",
    "- `return` is a keyword indicating the return value, as is often unnecessary  \n",
    "\n",
    "\n",
    "Let us make this example slightly better by “remembering” that `randn` can return a vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# still poor style\n",
    "function generatedata(n)\n",
    "    ϵ = randn(n) # use built in function\n",
    "\n",
    "    for i in eachindex(ϵ)\n",
    "        ϵ[i] = ϵ[i]^2 # squaring the result\n",
    "    end\n",
    "\n",
    "    return ϵ\n",
    "end\n",
    "data = generatedata(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While better, the looping over the `i` index to square the results is difficult to read\n",
    "\n",
    "Instead of looping, we can **broadcast** the `^2` square function over a vector using a `.`\n",
    "\n",
    "To be clear, unlike Python, R, and MATLAB (to a lesser extent), the reason to drop the `for` is **not** for performance reasons, but rather because of code clarity\n",
    "\n",
    "Loops of this sort are at least as efficient as vectorized approach in compiled languages like Julia, so use a for loop if you think it makes the code more clear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# better style\n",
    "function generatedata(n)\n",
    "    ϵ = randn(n) # use built in function\n",
    "    return ϵ.^2\n",
    " end\n",
    "data = generatedata(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even drop the `function` if we define it on a single line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# good style\n",
    "generatedata(n) = randn(n).^2\n",
    "data = generatedata(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can broadcast any function, where squaring is only a special case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# good style\n",
    "f(x) = x^2 # simple square function\n",
    "generatedata(n) = f.(randn(n)) # uses broadcast for some function `f`\n",
    "data = generatedata(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a final – abstract – approach, we can make the `generatedata` function able to generically apply to a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "generatedata(n, gen) = gen.(randn(n)) # uses broadcast for some function `gen`\n",
    "\n",
    "f(x) = x^2 # simple square function\n",
    "data = generatedata(5, f) # applies f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whether this example is better or worse than the previous version depends on how it is used\n",
    "\n",
    "High degrees of abstraction and generality, e.g. passing in a function `f` in this case, can make code either clearer or more confusing, but Julia enables you to use these techniques **with no performance overhead**\n",
    "\n",
    "For this particular case, the clearest and most general solution is probably the simplest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# direct solution with broadcasting, and small user-defined function\n",
    "n = 100\n",
    "f(x) = x^2\n",
    "\n",
    "x = randn(n)\n",
    "plot(f.(x), label=\"x^2\")\n",
    "plot!(x, label=\"x\") # layer on the same plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While broadcasting above superficially looks like vectorizing functions in MATLAB, or Python ufuncs, it is much richer and built on core foundations of the language\n",
    "\n",
    "The other additional function `plot!` adds a graph to the existing plot\n",
    "\n",
    "This follows a general convention in Julia, where a function that modifies the arguments or a global state has a `!` at the end of its name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A Slightly More Useful Function\n",
    "\n",
    "Let’s make a slightly more useful function\n",
    "\n",
    "This function will be passed in a choice of probability distribution and respond by plotting a histogram of observations\n",
    "\n",
    "In doing so we’ll make use of the `Distributions` package, which we assume was instantiated above with the project\n",
    "\n",
    "Here’s the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using Distributions\n",
    "\n",
    "function plothistogram(distribution, n)\n",
    "    ϵ = rand(distribution, n)  # n draws from distribution\n",
    "    histogram(ϵ)\n",
    "end\n",
    "\n",
    "lp = Laplace()\n",
    "plothistogram(lp, 500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s have a casual discussion of how all this works while leaving technical details for later in the lectures\n",
    "\n",
    "First, `lp = Laplace()` creates an instance of a data type defined\n",
    "in the `Distributions` module that represents the Laplace distribution\n",
    "\n",
    "The name `lp` is bound to this value\n",
    "\n",
    "When we make the function call `plothistogram(lp, 500)` the code in the body\n",
    "of the function `plothistogram` is run with\n",
    "\n",
    "- the name `distribution` bound to the same value as `lp`  \n",
    "- the name `n` bound to the integer `500`  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A Mystery\n",
    "\n",
    "Now consider the function call `rand(distribution, n)`\n",
    "\n",
    "This looks like something of a mystery\n",
    "\n",
    "The function `rand()` is defined in the base library such that `rand(n)` returns `n` uniform random variables on $ [0, 1) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "rand(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, `distribution` points to a data type representing the Laplace distribution that has been defined in a third party package\n",
    "\n",
    "So how can it be that `rand()` is able to take this kind of value as an\n",
    "argument and return the output that we want?\n",
    "\n",
    "The answer in a nutshell is **multiple dispatch**, which Julia uses to implement **generic programming**\n",
    "\n",
    "This refers to the idea that functions in Julia can have different behavior\n",
    "depending on the particular arguments that they’re passed\n",
    "\n",
    "Hence in Julia we can take an existing function and give it a new behavior by defining how it acts on a new type of value\n",
    "\n",
    "The compiler knows which function definition to apply to in a given setting by looking at the types of the values the function is called on\n",
    "\n",
    "In Julia these alternative versions of a function are called **methods**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Variations on Fixed Points\n",
    "\n",
    "Take a mapping $ f : X \\to X $ for some set $ X $\n",
    "\n",
    "If there exists an $ x^* \\in X $ such that $ f(x^*) = x^* $, then $ x^* $: is called a “fixed point” of $ f $\n",
    "\n",
    "For our second example, we will start with a simple example of determining fixed points of a function\n",
    "\n",
    "The goal is to start with code in a MATLAB style, and move towards a more **Julian** style with high mathematical clarity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed Point Maps\n",
    "\n",
    "Consider the simple equation, where the scalars $ p,\\beta $ are given, and  $ v $ is the scalar we wish to solve for\n",
    "\n",
    "$$\n",
    "v = p + \\beta v\n",
    "$$\n",
    "\n",
    "Of course, in this simple example, with parameter restrictions this can be solved as $ v = p/(1 - \\beta) $\n",
    "\n",
    "Rearrange the equation in terms of a map $ f(x) : \\mathbb R \\to \\mathbb R $\n",
    "\n",
    "\n",
    "<a id='equation-fixed-point-map'></a>\n",
    "$$\n",
    "v = f(v) \\tag{1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "f(v) := p + \\beta v\n",
    "$$\n",
    "\n",
    "Therefore, a fixed point $ v^* $ of $ f(\\cdot) $ is a solution to the above problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### While Loops\n",
    "\n",
    "One approach to finding a fixed point of [(1)](#equation-fixed-point-map) is to start with an initial value, and iterate the map\n",
    "\n",
    "\n",
    "<a id='equation-fixed-point-naive'></a>\n",
    "$$\n",
    "v^{n+1} = f(v^n) \\tag{2}\n",
    "$$\n",
    "\n",
    "For this exact `f` function,  we can see the convergence to $ v = p/(1-\\beta) $ when $ |\\beta| < 1 $ by iterating backwards and taking $ n\\to\\infty $\n",
    "\n",
    "$$\n",
    "v^{n+1} = p + \\beta v^n = p + \\beta p + \\beta^2 v^{n-1} = p \\sum_{i=0}^{n-1} \\beta^i + \\beta^n v_0\n",
    "$$\n",
    "\n",
    "To implement the iteration in [(2)](#equation-fixed-point-naive), we start by solving this problem with a `while` loop\n",
    "\n",
    "The syntax for the while loop contains no surprises, and looks nearly identical to a MATLAB implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# poor style\n",
    "p = 1.0 # note 1.0 rather than 1\n",
    "β = 0.9\n",
    "maxiter = 1000\n",
    "tolerance = 1.0E-7\n",
    "v_iv = 0.8 # initial condition\n",
    "\n",
    "# setup the algorithm\n",
    "v_old = v_iv\n",
    "normdiff = Inf\n",
    "iter = 1\n",
    "while normdiff > tolerance && iter <= maxiter\n",
    "    v_new = p + β * v_old # the f(v) map\n",
    "    normdiff = norm(v_new - v_old)\n",
    "\n",
    "    # replace and continue\n",
    "    v_old = v_new\n",
    "    iter = iter + 1\n",
    "end\n",
    "println(\"Fixed point = $v_old, and |f(x) - x| = $normdiff in $iter iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `while` loop, like the `for` loop should only be used directly in Jupyter or the inside of a function\n",
    "\n",
    "Here, we have used the `norm` function (from the `LinearAlgebra` base library) to compare the values\n",
    "\n",
    "The other new function is the `println` with the string interpolation, which splices the value of an expression or variable prefixed by `\\$` into a string\n",
    "\n",
    "An alternative approach is to use a `for` loop, and check for convergence in each iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# setup the algorithm\n",
    "v_old = v_iv\n",
    "normdiff = Inf\n",
    "iter = 1\n",
    "for i in 1:maxiter\n",
    "    v_new = p + β * v_old # the f(v) map\n",
    "    normdiff = norm(v_new - v_old)\n",
    "    if normdiff < tolerance # check convergence\n",
    "        iter = i\n",
    "        break # converged, exit loop\n",
    "    end\n",
    "    # replace and continue\n",
    "    v_old = v_new\n",
    "end\n",
    "println(\"Fixed point = $v_old, and |f(x) - x| = $normdiff in $iter iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new feature there is `break` , which leaves a `for` or `while` loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using a Function\n",
    "\n",
    "The first problem with this setup is that it depends on being sequentially run – which can be easily remedied with a function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# better, but still poor style\n",
    "function v_fp(β, ρ, v_iv, tolerance, maxiter)\n",
    "    # setup the algorithm\n",
    "    v_old = v_iv\n",
    "    normdiff = Inf\n",
    "    iter = 1\n",
    "    while normdiff > tolerance && iter <= maxiter\n",
    "        v_new = p + β * v_old # the f(v) map\n",
    "        normdiff = norm(v_new - v_old)\n",
    "\n",
    "        # replace and continue\n",
    "        v_old = v_new\n",
    "        iter = iter + 1\n",
    "    end\n",
    "    return (v_old, normdiff, iter) # returns a tuple\n",
    "end\n",
    "\n",
    "# some values\n",
    "p = 1.0 # note 1.0 rather than 1\n",
    "β = 0.9\n",
    "maxiter = 1000\n",
    "tolerance = 1.0E-7\n",
    "v_initial = 0.8 # initial condition\n",
    "\n",
    "v_star, normdiff, iter = v_fp(β, p, v_initial, tolerance, maxiter)\n",
    "println(\"Fixed point = $v_star, and |f(x) - x| = $normdiff in $iter iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While better, there could still be improvements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Passing a Function\n",
    "\n",
    "The chief issue is that the algorithm (finding a fixed point) is reusable and generic, while the function we calculate `p + β * v` is specific to our problem\n",
    "\n",
    "A key feature of languages like Julia, is the ability to efficiently handle functions passed to other functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# better style\n",
    "function fixedpointmap(f, iv, tolerance, maxiter)\n",
    "    # setup the algorithm\n",
    "    x_old = iv\n",
    "    normdiff = Inf\n",
    "    iter = 1\n",
    "    while normdiff > tolerance && iter <= maxiter\n",
    "        x_new = f(x_old) # use the passed in map\n",
    "        normdiff = norm(x_new - x_old)\n",
    "        x_old = x_new\n",
    "        iter = iter + 1\n",
    "    end\n",
    "    return (x_old, normdiff, iter)\n",
    "end\n",
    "\n",
    "# define a map and parameters\n",
    "p = 1.0\n",
    "β = 0.9\n",
    "f(v) = p + β * v # note that p and β are used in the function!\n",
    "\n",
    "maxiter = 1000\n",
    "tolerance = 1.0E-7\n",
    "v_initial = 0.8 # initial condition\n",
    "\n",
    "v_star, normdiff, iter = fixedpointmap(f, v_initial, tolerance, maxiter)\n",
    "println(\"Fixed point = $v_star, and |f(x) - x| = $normdiff in $iter iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much closer, but there are still hidden bugs if the user orders the settings or returns types wrong"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Named Arguments and Return Values\n",
    "\n",
    "To enable this, Julia has two features:  named function parameters, and named tuples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# good style\n",
    "function fixedpointmap(f; iv, tolerance=1E-7, maxiter=1000)\n",
    "    # setup the algorithm\n",
    "    x_old = iv\n",
    "    normdiff = Inf\n",
    "    iter = 1\n",
    "    while normdiff > tolerance && iter <= maxiter\n",
    "        x_new = f(x_old) # use the passed in map\n",
    "        normdiff = norm(x_new - x_old)\n",
    "        x_old = x_new\n",
    "        iter = iter + 1\n",
    "    end\n",
    "    return (value = x_old, normdiff=normdiff, iter=iter) # A named tuple\n",
    "end\n",
    "\n",
    "# define a map and parameters\n",
    "p = 1.0\n",
    "β = 0.9\n",
    "f(v) = p + β * v # note that p and β are used in the function!\n",
    "\n",
    "sol = fixedpointmap(f, iv=0.8, tolerance=1.0E-8) # don't need to pass\n",
    "println(\"Fixed point = $(sol.value), and |f(x) - x| = $(sol.normdiff) in $(sol.iter)\"*\n",
    "        \" iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, all function parameters after the `;` in the list, must be called by name\n",
    "\n",
    "Furthermore, a default value may be enabled – so the named parameter `iv` is required while `tolerance` and `maxiter` have default values\n",
    "\n",
    "The return type of the function also has named fields, `value, normdiff,` and `iter` – all accessed intuitively using `.`\n",
    "\n",
    "To show the flexibilty of this code, we can use it to find a fixed point of the non-linear logistic equation, $ x = f(x) $ where $ f(x) := r x (1-x) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "r = 2.0\n",
    "f(x) = r * x * (1 - x)\n",
    "\n",
    "sol = fixedpointmap(f, iv=0.8)\n",
    "println(\"Fixed point = $(sol.value), and |f(x) - x| = $(sol.normdiff) in $(sol.iter) iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using a Package\n",
    "\n",
    "But best of all is to avoid writing code altogether"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# best style\n",
    "using NLsolve\n",
    "\n",
    "p = 1.0\n",
    "β = 0.9\n",
    "f(v) = p .+ β * v # broadcast the +\n",
    "sol = fixedpoint(f, [0.8])\n",
    "println(\"Fixed point = $(sol.zero), and |f(x) - x| = $(norm(f(sol.zero) - sol.zero)) in \" *\n",
    "        \"$(sol.iterations) iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `fixedpoint` function from the `NLsolve.jl` library implements the simple fixed point iteration scheme above\n",
    "\n",
    "Since the `NLsolve` library only accepts vector based inputs, we needed to make the `f(v)` function broadcast on the `+` sign, and pass in the initial condition as a vector of length 1 with `[0.8]`\n",
    "\n",
    "While a key benefit of using a package is that the code is clearer, and the implementation is tested, by using an orthogonal library we also enable performance improvements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# best style\n",
    "p = 1.0\n",
    "β = 0.9\n",
    "iv = [0.8]\n",
    "sol = fixedpoint(v -> p .+ β * v, iv)\n",
    "println(\"Fixed point = $(sol.zero), and |f(x) - x| = $(norm(f(sol.zero) - sol.zero)) in \" *\n",
    "        \"$(sol.iterations) iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this completes in `3` iterations vs `177` for the naive fixed point iteration algorithm\n",
    "\n",
    "Since Anderson iteration is doing more calculations in an iteration,  whether it is faster or not would depend on the complexity of the `f` function\n",
    "\n",
    "But this demonstrates the value of keeping the math separate from the algorithm, since by decoupling the mathematical definition of the fixed point from the implementation in [(2)](#equation-fixed-point-naive), we were able to exploit new algorithms for finding a fixed point\n",
    "\n",
    "The only other change in this function is the move from directly defining `f(v)` and using an **anonymous** function\n",
    "\n",
    "Similar to anonymous functions in MATLAB, and lambda functions in Python, Julia enables the creation of small functions without any names\n",
    "\n",
    "The code `v -> p .+ β * v` defines a function of a dummy argument, `v` with the same body as our `f(x)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Composing Packages\n",
    "\n",
    "A key benefit of using Julia is that you can compose various packages, types, and techniques, without making changes to your underlying source\n",
    "\n",
    "As an example, consider if we want to solve the model with a higher-precision, as floating points cannot be distinguished beyond the machine epsilon for that type (recall that computers approximate real numbers to the nearest binary  of a given precision; the *machine epsilon* is the smallest nonzero magnitude)\n",
    "\n",
    "In Julia, this number can be calculated as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eps()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For many cases, this is sufficient precision – but consider that in iterative algorithms applied millions of times, those small differences can add up\n",
    "\n",
    "The only change we will need to our model in order to use a different floating point type is to call the function with an arbitrary precision floating point, `BigFloat`, for the initial value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# use arbitrary precision floating points\n",
    "p = 1.0\n",
    "β = 0.9\n",
    "iv = [BigFloat(0.8)] # higher precision\n",
    "\n",
    "# otherwise identical\n",
    "sol = fixedpoint(v -> p .+ β * v, iv)\n",
    "println(\"Fixed point = $(sol.zero), and |f(x) - x| = $(norm(f(sol.zero) - sol.zero)) in \" *\n",
    "        \"$(sol.iterations) iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the literal `BigFloat(0.8)` takes the number `0.8` and changes it to an arbitrary precision number\n",
    "\n",
    "The result is that the residual is now **exactly** `0.0` since it is able to use arbitrary precision in the calculations, and the solution has a finite-precision solution with those parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multivariate Fixed Point Maps\n",
    "\n",
    "The above example can be extended to multivariate maps without any modifications to the fixed point iteration code\n",
    "\n",
    "Using our own, homegrown iteration and simply passing in a bivariate map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p = [1.0, 2.0]\n",
    "β = 0.9\n",
    "iv = [0.8, 2.0]\n",
    "f(v) = p .+ β * v # note that p and β are used in the function!\n",
    "\n",
    "sol = fixedpointmap(f, iv = iv, tolerance = 1.0E-8)\n",
    "println(\"Fixed point = $(sol.value), and |f(x) - x| = $(sol.normdiff) in $(sol.iter)\"*\n",
    "\"iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also works without any modifications with the `fixedpoint` library function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using NLsolve\n",
    "\n",
    "p = [1.0, 2.0, 0.1]\n",
    "β = 0.9\n",
    "iv =[0.8, 2.0, 51.0]\n",
    "f(v) = p .+ β * v\n",
    "\n",
    "sol = fixedpoint(v -> p .+ β * v, iv)\n",
    "println(\"Fixed point = $(sol.zero), and |f(x) - x| = $(norm(f(sol.zero) - sol.zero)) in \" *\n",
    "        \"$(sol.iterations) iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, to demonstrate the importance of composing different libraries, use a `StaticArrays.jl` type, which provides an efficient implementation for small arrays and matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using NLsolve, StaticArrays\n",
    "p = @SVector [1.0, 2.0, 0.1]\n",
    "β = 0.9\n",
    "iv = @SVector  [0.8, 2.0, 51.0]\n",
    "f(v) = p .+ β * v\n",
    "\n",
    "sol = fixedpoint(v -> p .+ β * v, iv)\n",
    "println(\"Fixed point = $(sol.zero), and |f(x) - x| = $(norm(f(sol.zero) - sol.zero)) in \" *\n",
    "        \"$(sol.iterations) iterations\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `@SVector` in front of the `[1.0, 2.0, 0.1]` is a macro for turning a vector literal into a static vector\n",
    "\n",
    "All macros in Julia are prefixed by `@` in the name, and manipulate the code prior to compilation\n",
    "\n",
    "We will see a variety of macros, and discuss the “metaprogramming” behind them in a later lecture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='jbe-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Recall that $ n! $ is read as “$ n $ factorial” and defined as\n",
    "$ n! = n \\times (n - 1) \\times \\cdots \\times 2 \\times 1 $\n",
    "\n",
    "In Julia you can compute this value with `factorial(n)`\n",
    "\n",
    "Write your own version of this function, called `factorial2`, using a `for` loop\n",
    "\n",
    "\n",
    "<a id='jbe-ex2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "The [binomial random variable](https://en.wikipedia.org/wiki/Binomial_distribution) $ Y \\sim Bin(n, p) $ represents\n",
    "\n",
    "- number of successes in $ n $ binary trials  \n",
    "- each trial succeeds with probability $ p $  \n",
    "\n",
    "\n",
    "Using only `rand()` from the set of Julia’s built-in random number\n",
    "generators (not the `Distributions` package), write a function `binomial_rv` such that `binomial_rv(n, p)` generates one draw of $ Y $\n",
    "\n",
    "Hint: If $ U $ is uniform on $ (0, 1) $ and $ p \\in (0,1) $, then the expression `U < p` evaluates to `true` with probability $ p $\n",
    "\n",
    "\n",
    "<a id='jbe-ex3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    "Compute an approximation to $ \\pi $ using Monte Carlo\n",
    "\n",
    "For random number generation use only `rand()`\n",
    "\n",
    "Your hints are as follows:\n",
    "\n",
    "- If $ U $ is a bivariate uniform random variable on the unit square $ (0, 1)^2 $, then the probability that $ U $ lies in a subset $ B $ of $ (0,1)^2 $ is equal to the area of $ B $  \n",
    "- If $ U_1,\\ldots,U_n $ are iid copies of $ U $, then, as $ n $ gets larger, the fraction that falls in $ B $ converges to the probability of landing in $ B $  \n",
    "- For a circle, area = π * $ radius^2 $  \n",
    "\n",
    "\n",
    "\n",
    "<a id='jbe-ex4'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "\n",
    "Write a program that prints one realization of the following random device:\n",
    "\n",
    "- Flip an unbiased coin 10 times  \n",
    "- If 3 consecutive heads occur one or more times within this sequence, pay one dollar  \n",
    "- If not, pay nothing  \n",
    "\n",
    "\n",
    "Once again use only `rand()` as your random number generator\n",
    "\n",
    "\n",
    "<a id='jbe-ex5'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\n",
    "\n",
    "Simulate and plot the correlated time series\n",
    "\n",
    "$$\n",
    "x_{t+1} = \\alpha \\, x_t + \\epsilon_{t+1}\n",
    "\\quad \\text{where} \\quad\n",
    "x_0 = 0\n",
    "\\quad \\text{and} \\quad t = 0,\\ldots,n\n",
    "$$\n",
    "\n",
    "The sequence of shocks $ \\{\\epsilon_t\\} $ is assumed to be iid and standard normal\n",
    "\n",
    "Set $ n = 200 $ and $ \\alpha = 0.9 $\n",
    "\n",
    "\n",
    "<a id='jbe-ex6'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6\n",
    "\n",
    "Plot three simulated time series, one for each of the cases $ \\alpha = 0 $, $ \\alpha = 0.8 $ and $ \\alpha = 0.98 $\n",
    "\n",
    "(The figure will illustrate how time series with the same one-step-ahead conditional volatilities, as these three processes have, can have very different unconditional volatilities)\n",
    "\n",
    "\n",
    "<a id='jbe-ex7'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7\n",
    "\n",
    "This exercise is more challenging\n",
    "\n",
    "Take a random walk, starting from $ x_0 = 1 $\n",
    "\n",
    "$$\n",
    "x_{t+1} = \\, \\alpha \\, x_t + \\sigma\\, \\epsilon_{t+1}\n",
    "\\quad \\text{where} \\quad\n",
    "x_0 = 1\n",
    "\\quad \\text{and} \\quad t = 0,\\ldots,t_{\\max}\n",
    "$$\n",
    "\n",
    "- Furthermore, assume that the $ x_{t_{\\max}} = 0 $  (i.e. at $ t_{\\max} $, the value drops to zero, regardless of its current state)  \n",
    "- The sequence of shocks $ \\{\\epsilon_t\\} $ is assumed to be iid and standard normal  \n",
    "- For a given path $ \\{x_t\\} $ define a **first-passage time** as $ T_a = \\min\\{t\\, |\\, x_t \\leq a\\} $, where by the assumption of the process $ T_a \\leq t_{\\max} $  \n",
    "\n",
    "\n",
    "Start with $ \\sigma = 0.2, \\alpha = 1.0 $\n",
    "\n",
    "1. calculate the first-passage time, $ T_0 $, for 100 simulated random walks – to a $ t_{\\max} = 200 $ and plot a histogram  \n",
    "1. plot the sample mean of $ T_0 $ from the simulation for $ \\alpha \\in \\{0.8, 1.0, 1.2\\} $  \n",
    "\n",
    "\n",
    "\n",
    "<a id='jbe-ex8a'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 8(a)\n",
    "\n",
    "This exercise is more challenging\n",
    "\n",
    "The root of a univariate function $ f(\\cdot) $ is an $ x $ such that $ f(x) = 0 $\n",
    "\n",
    "One solution method to find local roots of smooth functions is called Newton’s method\n",
    "\n",
    "Starting with an $ x_0 $ guess, a function $ f(\\cdot) $ and the first-derivative $ f'(\\cdot) $, the algorithm is to repeat\n",
    "\n",
    "$$\n",
    "x^{n+1} = x^n - \\frac{f(x^n)}{f'(x^n)}\n",
    "$$\n",
    "\n",
    "until $ | x^{n+1} - x^n| $ is below a tolerance\n",
    "\n",
    "1. Use a variation of the `fixedpointmap` code to implement Newton’s method, where the function would accept arguments `f, f_prime, x_0, tolerance, maxiter`  \n",
    "1. Test it with $ f(x) = (x-1)^3 $ and another function of your choice where you can analytically find the derivative  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 8(b)\n",
    "\n",
    "For those impatient to use more advanced features of Julia, implement a version of Exercise 8(a) where `f_prime` is calculated with auto-differentiation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using ForwardDiff\n",
    "\n",
    "# operator to get the derivative of this function using AD\n",
    "D(f) = x -> ForwardDiff.derivative(f, x)\n",
    "\n",
    "# example usage: create a function and get the derivative\n",
    "f(x) = x^2\n",
    "f_prime = D(f)\n",
    "\n",
    "f(0.1), f_prime(0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Using the `D(f)` operator definition above, implement a version of Newton’s method that does not require the user to provide an analytical derivative  \n",
    "1. Test the sorts of `f` functions which can be automatically integrated by `ForwardDff.jl`  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "function factorial2(n)\n",
    "    k = 1\n",
    "    for i in 1:n\n",
    "        k *= i  # or k = k * i\n",
    "    end\n",
    "    return k\n",
    "end\n",
    "\n",
    "factorial2(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "factorial2(4) == factorial(4) # built-in function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "function binomial_rv(n, p)\n",
    "    count = 0\n",
    "    U = rand(n)\n",
    "    for i in 1:n\n",
    "        if U[i] < p\n",
    "            count += 1 # or count = count + 1\n",
    "        end\n",
    "    end\n",
    "    return count\n",
    "end\n",
    "\n",
    "for j in 1:25\n",
    "    b = binomial_rv(10, 0.5)\n",
    "    print(\"$b, \")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    "Consider a circle with diameter 1 embedded in a unit square\n",
    "\n",
    "Let $ A $ be its area and let $ r = 1/2 $ be its radius\n",
    "\n",
    "If we know $ \\pi $ then we can compute $ A $ via\n",
    "$ A = \\pi r^2 $\n",
    "\n",
    "But the point here is to compute $ \\pi $, which we can do by\n",
    "$ \\pi = A / r^2 $\n",
    "\n",
    "Summary: If we can estimate the area of the unit circle, then dividing\n",
    "by $ r^2 = (1/2)^2 = 1/4 $ gives an estimate of $ \\pi $\n",
    "\n",
    "We estimate the area by sampling bivariate uniforms and looking at the\n",
    "fraction that fall into the unit circle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = 1000000\n",
    "count = 0\n",
    "for i in 1:n\n",
    "    u, v = rand(2)\n",
    "    d = sqrt((u - 0.5)^2 + (v - 0.5)^2)  # distance from middle of square\n",
    "    if d < 0.5\n",
    "        count += 1\n",
    "    end\n",
    "end\n",
    "\n",
    "area_estimate = count / n\n",
    "\n",
    "print(area_estimate * 4)  # dividing by radius**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "payoff = 0\n",
    "count = 0\n",
    "\n",
    "print(\"Count = \")\n",
    "\n",
    "for i in 1:10\n",
    "    U = rand()\n",
    "    if U < 0.5\n",
    "        count += 1\n",
    "    else\n",
    "        count = 0\n",
    "    end\n",
    "    print(count)\n",
    "    if count == 3\n",
    "        payoff = 1\n",
    "    end\n",
    "end\n",
    "println(\"\\npayoff = $payoff\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can simplify this somewhat using the **ternary operator**. Here are\n",
    "some examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "a = 1  < 2 ? \"foo\" : \"bar\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "a = 1 > 2 ? \"foo\" : \"bar\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this construction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "payoff = 0.0\n",
    "count = 0.0\n",
    "\n",
    "print(\"Count = \")\n",
    "\n",
    "for i in 1:10\n",
    "    U = rand()\n",
    "    count = U < 0.5 ? count + 1 : 0\n",
    "    print(count)\n",
    "    if count == 3\n",
    "        payoff = 1\n",
    "    end\n",
    "end\n",
    "println(\"\\npayoff = $payoff\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\n",
    "\n",
    "Here’s one solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using Plots\n",
    "gr(fmt=:png); # setting for easier display in jupyter notebooks\n",
    "α = 0.9\n",
    "n = 200\n",
    "x = zeros(n + 1)\n",
    "\n",
    "for t in 1:n\n",
    "    x[t+1] = α * x[t] + randn()\n",
    "end\n",
    "plot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "αs = [0.0, 0.8, 0.98]\n",
    "n = 200\n",
    "p = plot() # naming a plot to add to\n",
    "\n",
    "for α in αs\n",
    "    x = zeros(n + 1)\n",
    "    x[1] = 0.0\n",
    "    for t in 1:n\n",
    "        x[t+1] = α * x[t] + randn()\n",
    "    end\n",
    "    plot!(p, x, label = \"alpha = $α\") # add to plot p\n",
    "end\n",
    "p # display plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7: Hint\n",
    "\n",
    "As a hint, notice the following pattern for finding the number of draws of a uniform random number until it is below a given threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "function drawsuntilthreshold(threshold; maxdraws=100)\n",
    "    for i in 1:maxdraws\n",
    "        val = rand()\n",
    "        if val < threshold # checks threshold\n",
    "            return i # leaves function, returning draw number\n",
    "        end\n",
    "    end\n",
    "    return Inf # if here, reached maxdraws\n",
    "end\n",
    "\n",
    "draws = drawsuntilthreshold(0.2, maxdraws=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, it is sometimes convenient to add to just push numbers onto an array without indexing it directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "vals = zeros(0) # empty vector\n",
    "\n",
    "for i in 1:100\n",
    "    val = rand()\n",
    "    if val < 0.5\n",
    "        push!(vals, val)\n",
    "    end\n",
    "end\n",
    "println(\"There were $(length(vals)) below 0.5\")"
   ]
  }
 ],
 "metadata": {
  "filename": "julia_by_example.rst",
  "kernelspec": {
   "display_name": "Julia 1.2",
   "language": "julia",
   "name": "julia-1.2"
  },
  "title": "Introductory Examples"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}