{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing integral with quasi-Monte Carlo methods\n",
    "\n",
    "**Randall Romero Aguilar, PhD**\n",
    "\n",
    "This demo is based on the original Matlab demo accompanying the  <a href=\"https://mitpress.mit.edu/books/applied-computational-economics-and-finance\">Computational Economics and Finance</a> 2001 textbook by Mario Miranda and Paul Fackler.\n",
    "\n",
    "Original (Matlab) CompEcon file: **demqua01bis.m**\n",
    "\n",
    "Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
    "\n",
    "    !pip install compecon --upgrade\n",
    "\n",
    "<i>Last updated: 2021-Oct-01</i>\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## About"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To seven significant digits,\n",
    "\\begin{align*}\n",
    "A &=\\int_{-1}^1\\int_{-1}^1  e^{-x_1}\\cos^2(x_2)dx _1dx_2\\\\\n",
    " &=\\int_{-1}^1 e^{-x_1} dx _1 \\times \\int_{-1}^1 \\cos^2(x_2) dx_2\\\\\n",
    " &=\\left(e - \\tfrac{1}{e}\\right) \\times \\left(1+\\tfrac{1}{2}\\sin(2)\\right)\\\\\n",
    " &\\approx 3.4190098\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial tasks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from compecon import qnwequi\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make support function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = lambda x1: np.exp(-x1)\n",
    "f2 = lambda x2: np.cos(x2)**2\n",
    "f = lambda x1, x2: f1(x1) * f2(x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quad(method, n):\n",
    "    (x1, x2), w = qnwequi(n,[-1, -1], [1, 1],method)\n",
    "    return w.dot(f(x1, x2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the approximation errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "nlist = range(3,7)\n",
    "quadmethods = ['Random', 'Neiderreiter','Weyl']\n",
    "\n",
    "f_quad = np.array([[quad(qnw[0], 10**ni) for qnw in quadmethods] for ni in nlist])\n",
    "f_true = (np.exp(1) - np.exp(-1)) * (1+0.5*np.sin(2))\n",
    "f_error = np.log10(np.abs(f_quad/f_true - 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make table with results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": "          Random  Neiderreiter      Weyl\nNodes                                   \n$10^3$ -1.804202     -2.879952 -3.518213\n$10^4$ -2.603463     -3.148008 -3.929531\n$10^5$ -2.837374     -3.982252 -4.377373\n$10^6$ -3.617125     -5.449389 -5.740493",
      "text/html": "<div>\n<style scoped>\n    .dataframe tbody tr th:only-of-type {\n        vertical-align: middle;\n    }\n\n    .dataframe tbody tr th {\n        vertical-align: top;\n    }\n\n    .dataframe thead th {\n        text-align: right;\n    }\n</style>\n<table border=\"1\" class=\"dataframe\">\n  <thead>\n    <tr style=\"text-align: right;\">\n      <th></th>\n      <th>Random</th>\n      <th>Neiderreiter</th>\n      <th>Weyl</th>\n    </tr>\n    <tr>\n      <th>Nodes</th>\n      <th></th>\n      <th></th>\n      <th></th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>$10^3$</th>\n      <td>-1.804202</td>\n      <td>-2.879952</td>\n      <td>-3.518213</td>\n    </tr>\n    <tr>\n      <th>$10^4$</th>\n      <td>-2.603463</td>\n      <td>-3.148008</td>\n      <td>-3.929531</td>\n    </tr>\n    <tr>\n      <th>$10^5$</th>\n      <td>-2.837374</td>\n      <td>-3.982252</td>\n      <td>-4.377373</td>\n    </tr>\n    <tr>\n      <th>$10^6$</th>\n      <td>-3.617125</td>\n      <td>-5.449389</td>\n      <td>-5.740493</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results = pd.DataFrame(f_error, columns=quadmethods)\n",
    "results['Nodes'] = ['$10^%d$' % n for n in nlist]\n",
    "results.set_index('Nodes', inplace=True)\n",
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#results.to_latex('demqua01bis.tex', escape=False, float_format='%.1f')"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}