_notebooks/2020-06-22-02-Advanced-Applications-of-Simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Applications of Simulation\n",
    "> In this chapter, students will be introduced to some basic and advanced applications of simulation to solve real-world problems. We'll work through a business planning problem, learn about Monte Carlo Integration, Power Analysis with simulation and conclude with a financial portfolio simulation. After completing this chapter, students will be ready to apply simulation to solve everyday problems. This is the Summary of lecture \"Statistical Simulation in Python\", via datacamp.\n",
    "\n",
    "- toc: true \n",
    "- badges: true\n",
    "- comments: true\n",
    "- author: Chanseok Kang\n",
    "- categories: [Python, Datacamp, Statistics, Modeling]\n",
    "- image: images/power_analysis.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (10, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation for Business Planning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling Corn Production\n",
    "Suppose that you manage a small corn farm and are interested in optimizing your costs. In this exercise, we will model the production of corn.\n",
    "\n",
    "For simplicity, let's assume that corn production depends on only two factors: rain, which you don't control, and cost, which you control. Rain is normally distributed with mean 50 and standard deviation 15. For now, let's fix cost at 5,000. Corn produced in any season is a Poisson random variable while the average corn production is governed by the equation:\n",
    "\n",
    "$$ 100 * (cost)^{0.1} * (rain)^{0.2} $$\n",
    "Let's model this production function and simulate one outcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated Corn Production = 560\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(223)\n",
    "\n",
    "# Initialize variables\n",
    "cost = 5000\n",
    "rain = np.random.normal(50, 15)\n",
    "\n",
    "# Corn Production Model\n",
    "def corn_produced(rain, cost):\n",
    "    mean_corn = np.mean(100 * cost ** 0.1 * rain ** 0.2)\n",
    "    corn = np.random.poisson(mean_corn)\n",
    "    return corn\n",
    "\n",
    "# Simulate and print corn production\n",
    "corn_result = corn_produced(rain, cost)\n",
    "print(\"Simulated Corn Production = {}\".format(corn_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling Profits\n",
    "In the previous exercise, you built a model of corn production. For a small farm, you typically have no control over the price or demand for corn. Suppose that price is normally distributed with mean 40 and standard deviation 10. You are given a function `corn_demanded()`, which takes the price and determines the demand for corn. This is reasonable because demand is usually determined by the market and is not in your control.\n",
    "\n",
    "In this exercise, you will work on a function to calculate the profit by pulling together all the other simulated variables. The only input to this function will be the cost. Upon completion, you will have a function that will give you one simulated profit outcome for a given cost. This function can then be used for planning your costs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def corn_demanded(price):\n",
    "    mean_corn = 1000 - 8 * price\n",
    "    corn = np.random.poisson(abs(mean_corn))\n",
    "    return corn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated profit = 20675.3291075312\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(223)\n",
    "\n",
    "# Function to calculate profits\n",
    "def profits(cost):\n",
    "    rain = np.random.normal(50, 15)\n",
    "    price = np.random.normal(40, 10)\n",
    "    supply = corn_produced(rain, cost)\n",
    "    demand = corn_demanded(price)\n",
    "    equil_short = supply <= demand\n",
    "    if equil_short == True:\n",
    "        tmp = supply * price - cost\n",
    "        return tmp\n",
    "    else: \n",
    "        tmp2 = demand * price - cost\n",
    "        return tmp2\n",
    "result = profits(cost)\n",
    "print(\"Simulated profit = {}\".format(result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimizing Costs\n",
    "Now we will use the functions you've built to optimize our cost of production. We are interested in maximizing average profits. However, our profits depend on a number of factors, but we only control cost. Thus, we can simulate the uncertainty in the other factors and vary cost to see how our profits are impacted.\n",
    "\n",
    "Since you manage the small corn farm, you have the ability to choose your cost - from \\\\$100  to \\\\$5,000. You want to choose the cost that gives you the maximum average profit. In this exercise, we will simulate multiple outcomes for each cost level and calculate an average. We will then choose the cost that gives us the maximum mean profit. Upon completion, you will have a framework for selecting optimal inputs for business decisions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/chanseok/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average profit is maximized when cost = 1400\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(573)\n",
    "\n",
    "# Initialize results and cost_levels variables\n",
    "sims, results = 1000, {}\n",
    "cost_levels = np.arange(100, 5100, 100)\n",
    "\n",
    "# For each cost level, simulate profits and store mean profit\n",
    "for cost in cost_levels:\n",
    "    tmp_profits = []\n",
    "    for i in range(sims):\n",
    "        tmp_profits.append(profits(cost))\n",
    "    results[cost] = np.mean(tmp_profits)\n",
    "    \n",
    "# Get the cost that maximizes average profit\n",
    "cost_max = [x for x in results.keys() if results[x] == max(results.values())][0]\n",
    "print(\"Average profit is maximized when cost = {}\".format(cost_max))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo Integration\n",
    "- Calculate overall area\n",
    "- Randomly sample points in the area\n",
    "- Multiply the fraction of the points below the curve by overall area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integrating a Simple Function\n",
    "This is a simple exercise introducing the concept of Monte Carlo Integration.\n",
    "\n",
    "Here we will evaluate a simple integral $\\int_0^1 x e^x dx$. We know that the exact answer is 1, but simulation will give us an approximate solution, so we can expect an answer close to 1. As we saw in the video, it's a simple process. For a function of a single variable $f(x)$:\n",
    "\n",
    "1. Get the limits of the x-axis ($x_{min},x_{max}$) and y-axis ($max(f(x)),min(min(f(x)),0)$).\n",
    "2. Generate a number of uniformly distributed point in this box.\n",
    "3. Multiply the area of the box ($(max(f(x)−min(f(x))×(x_{max}−x_{min})$) by the fraction of points that lie below $f(x)$.\n",
    "\n",
    "Upon completion, you will have a framework for handling definite integrals using Monte Carlo Integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated answer = 0.7240166789450252, Actual Answer = 1\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "\n",
    "# Define the sim_integrate function\n",
    "def sim_integrate(func, xmin, xmax, sims):\n",
    "    x = np.random.uniform(xmin, xmax, sims)\n",
    "    y = np.random.uniform(min(min(func(x)), 0), max(func(x)), sims)\n",
    "    area = (max(y) - min(y)) * (xmax - xmin)\n",
    "    result = area * sum(abs(y) < abs(func(x)))/sims\n",
    "    return result\n",
    "    \n",
    "# Call the sim_integrate function and print results\n",
    "result = sim_integrate(func = lambda x: x * np.exp(x), xmin = 0, xmax = 1, sims = 50)\n",
    "print(\"Simulated answer = {}, Actual Answer = 1\".format(result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the value of pi\n",
    "Now we work through a classic example - estimating the value of $\\pi$.\n",
    "\n",
    "Imagine a square of side 2 with the origin $(0,0)$ as its center and the four corners having coordinates $(1,1),(1,−1),(−1,1),(−1,−1)$. The area of this square is $2×2=4$. Now imagine a circle of radius 1 with its center at the origin fitting perfectly inside this square. The area of the circle will be $\\pi × \\text{radius}^2=\\pi$.\n",
    "\n",
    "To estimate π, we randomly sample multiple points in this square & get the fraction of points inside the circle ($x^2+y^2 \\leq 1$). The area of the circle then is 4 times this fraction, which gives us our estimate of $\\pi$.\n",
    "\n",
    "After this exercise, you'll have a grasp of how to use simulation for computation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulated value of pi = 3.1468\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "\n",
    "# Initialize sims and circle_points\n",
    "sims, circle_points = 10000, 0 \n",
    "\n",
    "for i in range(sims):\n",
    "    # Generate the two coordinates of a point\n",
    "    point = np.random.uniform(-1, 1, size=2)\n",
    "    # if the point lies within the unit circle, increment counter\n",
    "    within_circle = point[0]**2 + point[1]**2 <= 1\n",
    "    if within_circle == True:\n",
    "        circle_points +=1\n",
    "        \n",
    "# Estimate pi as 4 times the avg number of points in the circle.\n",
    "pi_sim = np.mean(circle_points) * 4 / sims\n",
    "print(\"Simulated value of pi = {}\".format(pi_sim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation for Power Analysis\n",
    "- What is power?\n",
    "    - Power = $P(\\text{rejecting Null } \\vert \\text{ true alternative})$\n",
    "    - Probability of detecting an effect if it exists\n",
    "    - Depends on sample size, $\\alpha$ and effect size.\n",
    "    - Typically 80% power recommended for $\\alpha = 0.05$\n",
    "- Simulation for Power Analysis\n",
    "![power](image/power_analysis.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power Analysis - Part I\n",
    "Now we turn to power analysis. You typically want to ensure that any experiment or A/B test you run has at least 80% power. One way to ensure this is to calculate the sample size required to achieve 80% power.\n",
    "\n",
    "Suppose that you are in charge of a news media website and you are interested in increasing the amount of time users spend on your website. Currently, the time users spend on your website is normally distributed with a mean of 1 minute and a variance of 0.5 minutes. Suppose that you are introducing a feature that loads pages faster and want to know the sample size required to measure a 10% increase in time spent on the website.\n",
    "\n",
    "In this exercise, we will set up the framework to run one simulation, run a t-test, & calculate the p-value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-value: 0.5766409395002308, Statistically Significant? False\n"
     ]
    }
   ],
   "source": [
    "import scipy.stats as st\n",
    "\n",
    "np.random.seed(123)\n",
    "\n",
    "# Initialize effect_size, control_mean, control_sd\n",
    "effect_size, sample_size, control_mean, control_sd = 0.05, 50, 1, 0.5\n",
    "\n",
    "# Simulate control_time_spent and treatment_time_spent, assuming equal variance\n",
    "control_time_spent = np.random.normal(loc=control_mean, scale=control_sd, size=sample_size)\n",
    "treatment_time_spent = np.random.normal(loc=control_mean*(1+effect_size), \n",
    "                                        scale=control_sd, size=sample_size)\n",
    "\n",
    "# Run the t-test and get the p_value\n",
    "t_stat, p_value = st.ttest_ind(control_time_spent, treatment_time_spent)\n",
    "stat_sig = p_value < 0.05\n",
    "print(\"P-value: {}, Statistically Significant? {}\".format(p_value, stat_sig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power Analysis - Part II\n",
    "Previously, we simulated one instance of the experiment & generated a p-value. We will now use this framework to calculate statistical power. Power of an experiment is the experiment's ability to detect a difference between treatment & control if the difference really exists. It's good statistical hygiene to strive for 80% power.\n",
    "\n",
    "For our website, we want to know how many people need to visit each variant, such that we can detect a 10% increase in time spent with 80% power. For this, we start with a small sample (50), simulate multiple instances of this experiment & check power. If 80% power is reached, we stop. If not, we increase the sample size & try again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For 80% power, sample size required = 1580\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "sample_size = 50\n",
    "\n",
    "# Keep incrementing sample size by 10 till we reach required power\n",
    "while 1:\n",
    "    control_time_spent = np.random.normal(loc=control_mean, \n",
    "                                          scale=control_sd, \n",
    "                                          size=(sample_size, sims))\n",
    "    treatment_time_spent = np.random.normal(loc=control_mean * (1 + effect_size), \n",
    "                                            scale=control_sd, \n",
    "                                            size=(sample_size, sims))\n",
    "    t, p = st.ttest_ind(treatment_time_spent, control_time_spent)\n",
    "    \n",
    "    # Power is the fraction of times in the simulation when the p-value was less than 0.05\n",
    "    power = (p < 0.05).sum() / sims\n",
    "    if power >= 0.8: \n",
    "        break\n",
    "    else: \n",
    "        sample_size += 10\n",
    "print(\"For 80% power, sample size required = {}\".format(sample_size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications in Finance\n",
    "- Option & Instrument Pricing\n",
    "- Project Finance\n",
    "- Portfolio Evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Portfolio Simulation - Part I\n",
    "In the next few exercises, you will calculate the expected returns of a stock portfolio & characterize its uncertainty.\n",
    "\n",
    "Suppose you have invested \\\\$10,000 in your portfolio comprising of multiple stocks. You want to evaluate the portfolio's performance over 10 years. You can tweak your overall expected rate of return and volatility (standard deviation of the rate of return). Assume the rate of return follows a normal distribution.\n",
    "\n",
    "First, let's write a function that takes the principal (initial investment), number of years, expected rate of return and volatility as inputs and returns the portfolio's total value after 10 years.\n",
    "\n",
    "Upon completion of this exercise, you will have a function you can call to determine portfolio performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Portfolio return after 5 years = 1021.4013412039292\n"
     ]
    }
   ],
   "source": [
    "# rates is a Normal random variable and has size equal to number of years\n",
    "def portfolio_return(yrs, avg_return, volatility, principal):\n",
    "    np.random.seed(123)\n",
    "    rates = np.random.normal(loc=avg_return, scale=volatility, size=yrs)\n",
    "    # Calculate the return at the end of the period\n",
    "    end_return = principal\n",
    "    for x in rates:\n",
    "        end_return = end_return * (1 + x)\n",
    "    return end_return\n",
    "\n",
    "result = portfolio_return(yrs = 5, avg_return = 0.07, volatility = 0.15, principal = 1000)\n",
    "print(\"Portfolio return after 5 years = {}\".format(result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Portfolio Simulation - Part II\n",
    "Now we will use the simulation function you built to evaluate 10-year returns.\n",
    "\n",
    "Your stock-heavy portfolio has an initial investment of \\\\$10,000, an expected return of 7% and a volatility of 30%. You want to get a 95% confidence interval of what your investment will be worth in 10 years. We will simulate multiple samples of 10-year returns and calculate the confidence intervals on the distribution of returns.\n",
    "\n",
    "By the end of this exercise, you will have run a complete portfolio simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% CI of Returns: Lower = 3859.345207073691, Upper = 3859.345207073691\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "# Run 1,000 iterations and store the results\n",
    "sims, rets = 1000, []\n",
    "\n",
    "for i in range(sims):\n",
    "    rets.append(portfolio_return(yrs = 10, avg_return = 0.07, \n",
    "                                 volatility = 0.3, principal = 10000))\n",
    "\n",
    "# Calculate the 95% CI\n",
    "lower_ci = np.percentile(rets, 2.5)\n",
    "upper_ci = np.percentile(rets, 97.5)\n",
    "print(\"95% CI of Returns: Lower = {}, Upper = {}\".format(lower_ci, upper_ci))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Portfolio Simulation - Part III\n",
    "Previously, we ran a complete simulation to get a distribution for 10-year returns. Now we will use simulation for decision making.\n",
    "\n",
    "Let's go back to your stock-heavy portfolio with an expected return of 7% and a volatility of 30%. You have the choice of rebalancing your portfolio with some bonds such that the expected return is 4% & volatility is 10%. You have a principal of \\\\$10,000. You want to select a strategy based on how much your portfolio will be worth in 10 years. Let's simulate returns for both the portfolios and choose based on the least amount you can expect with 75% probability (25th percentile).\n",
    "\n",
    "Upon completion, you will know how to use a portfolio simulation for investment decisions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rets_stock = []\n",
    "rets_bond = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sticking to stocks gets you an additional return of -6696.359667655672\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(123)\n",
    "\n",
    "for i in range(sims):\n",
    "    rets_stock.append(portfolio_return(yrs = 10, avg_return = 0.07, \n",
    "                                       volatility = 0.3, principal = 10000))\n",
    "    rets_bond.append(portfolio_return(yrs = 10, avg_return = 0.04, \n",
    "                                      volatility = 0.1, principal = 10000))\n",
    "\n",
    "# Calculate the 25th percentile of the distributions and the amount you'd lose or gain\n",
    "rets_stock_perc = np.percentile(rets_stock, 25)\n",
    "rets_bond_perc = np.percentile(rets_bond, 25)\n",
    "additional_returns = rets_stock_perc - rets_bond_perc\n",
    "print(\"Sticking to stocks gets you an additional return of {}\".format(additional_returns))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}