Chapter 02 - The Importance of Linearity in Finance.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The capital asset pricing model and the security market line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "Linear regression with SciPy \n",
    "\"\"\"\n",
    "from scipy import stats\n",
    "\n",
    "stock_returns = [0.065, 0.0265, -0.0593, -0.001, 0.0346]\n",
    "mkt_returns = [0.055, -0.09, -0.041, 0.045, 0.022]\n",
    "beta, alpha, r_value, p_value, std_err = \\\n",
    "    stats.linregress(stock_returns, mkt_returns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5077431878770808 -0.008481900352462384\n"
     ]
    }
   ],
   "source": [
    "print(beta, alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multivariate linear regression of factor models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "Least squares regression with statsmodels \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "\n",
    "# Generate some sample data\n",
    "num_periods = 9\n",
    "all_values = np.array([np.random.random(8) \\\n",
    "                       for i in range(num_periods)])\n",
    "\n",
    "# Filter the data\n",
    "y_values = all_values[:, 0] # First column values as Y\n",
    "x_values = all_values[:, 1:] # All other values as X\n",
    "x_values = sm.add_constant(x_values) # Include the intercept\n",
    "results = sm.OLS(y_values, x_values).fit() # Regress and fit the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       1.000\n",
      "Model:                            OLS   Adj. R-squared:                  0.999\n",
      "Method:                 Least Squares   F-statistic:                     1117.\n",
      "Date:                Sun, 19 Aug 2018   Prob (F-statistic):             0.0230\n",
      "Time:                        00:30:24   Log-Likelihood:                 38.819\n",
      "No. Observations:                   9   AIC:                            -61.64\n",
      "Df Residuals:                       1   BIC:                            -60.06\n",
      "Df Model:                           7                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.1570      0.021      7.335      0.086      -0.115       0.429\n",
      "x1            -0.4513      0.014    -33.118      0.019      -0.624      -0.278\n",
      "x2             0.8468      0.012     68.712      0.009       0.690       1.003\n",
      "x3             0.1357      0.020      6.927      0.091      -0.113       0.385\n",
      "x4            -0.5728      0.020    -29.172      0.022      -0.822      -0.323\n",
      "x5             0.9883      0.025     39.539      0.016       0.671       1.306\n",
      "x6             0.1418      0.016      8.663      0.073      -0.066       0.350\n",
      "x7            -0.4438      0.015    -29.940      0.021      -0.632      -0.255\n",
      "==============================================================================\n",
      "Omnibus:                        0.482   Durbin-Watson:                   1.706\n",
      "Prob(Omnibus):                  0.786   Jarque-Bera (JB):                0.027\n",
      "Skew:                          -0.074   Prob(JB):                        0.987\n",
      "Kurtosis:                       2.778   Cond. No.                         19.2\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\python37\\lib\\site-packages\\scipy\\stats\\stats.py:1394: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=9\n",
      "  \"anyway, n=%i\" % int(n))\n"
     ]
    }
   ],
   "source": [
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.15700937 -0.45126648  0.84676324  0.13568198 -0.57281289  0.98829748\n",
      "  0.14180069 -0.44376925]\n"
     ]
    }
   ],
   "source": [
    "print(results.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A maximization example with linear programming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\" \n",
    "A linear optimization problem with 2 variables\n",
    "\"\"\"\n",
    "import pulp\n",
    "\n",
    "x = pulp.LpVariable('x', lowBound=0)\n",
    "y = pulp.LpVariable('y', lowBound=0)\n",
    "\n",
    "problem = pulp.LpProblem(\n",
    "    'A simple maximization objective', \n",
    "    pulp.LpMaximize)\n",
    "problem += 3*x + 2*y, 'The objective function'\n",
    "problem += 2*x + y <= 100, '1st constraint'\n",
    "problem += x + y <= 80, '2nd constraint'\n",
    "problem += x <= 40, '3rd constraint'\n",
    "problem.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximization Results:\n",
      "x = 20.0\n",
      "y = 60.0\n"
     ]
    }
   ],
   "source": [
    "print(\"Maximization Results:\")\n",
    "for variable in problem.variables():\n",
    "    print(variable.name, '=', variable.varValue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A minimization example with integer programming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "An example of implementing an integer \n",
    "programming model with binary conditions \n",
    "\"\"\"\n",
    "import pulp\n",
    "\n",
    "dealers = ['X', 'Y', 'Z']\n",
    "variable_costs = {'X': 500, 'Y': 350, 'Z': 450}\n",
    "fixed_costs = {'X': 4000, 'Y': 2000, 'Z': 6000}\n",
    "\n",
    "# Define PuLP variables to solve\n",
    "quantities = pulp.LpVariable.dicts('quantity', \n",
    "                                   dealers, \n",
    "                                   lowBound=0,\n",
    "                                   cat=pulp.LpInteger)\n",
    "is_orders = pulp.LpVariable.dicts('orders', \n",
    "                                  dealers,\n",
    "                                  cat=pulp.LpBinary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Non-constant expressions cannot be multiplied",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-3-54295b53a785>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      6\u001b[0m model = pulp.LpProblem('A cost minimization problem',\n\u001b[0;32m      7\u001b[0m                        pulp.LpMinimize)\n\u001b[1;32m----> 8\u001b[1;33m \u001b[0mmodel\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvariable_costs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m*\u001b[0m                \u001b[0mquantities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m                \u001b[0mfixed_costs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mis_orders\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m               \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdealers\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'Minimize portfolio cost'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      9\u001b[0m \u001b[0mmodel\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mquantities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdealers\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m150\u001b[0m    \u001b[1;33m,\u001b[0m \u001b[1;34m'Total contracts required'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[0mmodel\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m30\u001b[0m \u001b[1;33m<=\u001b[0m \u001b[0mquantities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'X'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m<=\u001b[0m \u001b[1;36m100\u001b[0m    \u001b[1;33m,\u001b[0m \u001b[1;34m'Boundary of total volume of X'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-3-54295b53a785>\u001b[0m in \u001b[0;36m<listcomp>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m      6\u001b[0m model = pulp.LpProblem('A cost minimization problem',\n\u001b[0;32m      7\u001b[0m                        pulp.LpMinimize)\n\u001b[1;32m----> 8\u001b[1;33m \u001b[0mmodel\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvariable_costs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m*\u001b[0m                \u001b[0mquantities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m                \u001b[0mfixed_costs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mis_orders\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m               \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdealers\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'Minimize portfolio cost'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      9\u001b[0m \u001b[0mmodel\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mquantities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mdealers\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m150\u001b[0m    \u001b[1;33m,\u001b[0m \u001b[1;34m'Total contracts required'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m \u001b[0mmodel\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m30\u001b[0m \u001b[1;33m<=\u001b[0m \u001b[0mquantities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'X'\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m<=\u001b[0m \u001b[1;36m100\u001b[0m    \u001b[1;33m,\u001b[0m \u001b[1;34m'Boundary of total volume of X'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\python37\\lib\\site-packages\\pulp\\pulp.py\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m    785\u001b[0m                         \u001b[0me\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mc\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    786\u001b[0m         \u001b[1;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mLpVariable\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 787\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mLpAffineExpression\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    788\u001b[0m         \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    789\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mother\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mc:\\python37\\lib\\site-packages\\pulp\\pulp.py\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m    773\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    774\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 775\u001b[1;33m                     \u001b[1;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Non-constant expressions cannot be multiplied\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    776\u001b[0m                 \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    777\u001b[0m                     \u001b[0mc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconstant\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mTypeError\u001b[0m: Non-constant expressions cannot be multiplied"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "This is an example of implementing an integer programming model\n",
    "with binary variables the wrong way.\n",
    "\"\"\"\n",
    "# Initialize the model with constraints\n",
    "model = pulp.LpProblem('A cost minimization problem',\n",
    "                       pulp.LpMinimize)\n",
    "model += sum([(variable_costs[i] * \\\n",
    "               quantities[i] + \\\n",
    "               fixed_costs[i])*is_orders[i] \\\n",
    "              for i in dealers]), 'Minimize portfolio cost'\n",
    "model += sum([quantities[i] for i in dealers]) == 150\\\n",
    "    , 'Total contracts required'\n",
    "model += 30 <= quantities['X'] <= 100\\\n",
    "    , 'Boundary of total volume of X'\n",
    "model += 30 <= quantities['Y'] <= 90\\\n",
    "    , 'Boundary of total volume of Y'\n",
    "model += 30 <= quantities['Z'] <= 70\\\n",
    "    , 'Boundary of total volume of Z'\n",
    "model.solve()  # You will get an error running this code!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "This is an example of implementing an \n",
    "IP model with binary variables the correct way.\n",
    "\"\"\"\n",
    "# Initialize the model with constraints\n",
    "model = pulp.LpProblem('A cost minimization problem',\n",
    "                       pulp.LpMinimize)\n",
    "model += sum(\n",
    "    [variable_costs[i]*quantities[i] + \\\n",
    "         fixed_costs[i]*is_orders[i] for i in dealers])\\\n",
    "    , 'Minimize portfolio cost'\n",
    "model += sum([quantities[i] for i in dealers]) == 150\\\n",
    "    ,  'Total contracts required'\n",
    "model += is_orders['X']*30 <= quantities['X'] <= \\\n",
    "    is_orders['X']*100, 'Boundary of total volume of X'\n",
    "model += is_orders['Y']*30 <= quantities['Y'] <= \\\n",
    "    is_orders['Y']*90, 'Boundary of total volume of Y'\n",
    "model += is_orders['Z']*30 <= quantities['Z'] <= \\\n",
    "    is_orders['Z']*70, 'Boundary of total volume of Z'\n",
    "model.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimization Results:\n",
      "orders_X = 0.0\n",
      "orders_Y = 1.0\n",
      "orders_Z = 1.0\n",
      "quantity_X = 0.0\n",
      "quantity_Y = 90.0\n",
      "quantity_Z = 60.0\n",
      "Total cost: 66500.0\n"
     ]
    }
   ],
   "source": [
    "print('Minimization Results:')\n",
    "for variable in model.variables():\n",
    "    print(variable, '=', variable.varValue)\n",
    "\n",
    "print('Total cost:',  pulp.value(model.objective))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving linear equations using matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "Linear algebra with NumPy matrices \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "\n",
    "A = np.array([[2, 1, 1],[1, 3, 2],[1, 0, 0]])\n",
    "B = np.array([4, 5, 6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `linalg.solve` function of NumPy to solve a system of linear scalar\n",
    "equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  6.  15. -23.]\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.solve(A, B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The LU decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    " \"\"\" \n",
    "LU decomposition with SciPy \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import scipy.linalg as linalg\n",
    "\n",
    "\n",
    "# Define A and B\n",
    "A = np.array([\n",
    "    [2., 1., 1.],\n",
    "    [1., 3., 2.],\n",
    "    [1., 0., 0.]])\n",
    "B = np.array([4., 5., 6.])\n",
    "\n",
    "# Perform LU decompositiopn\n",
    "LU = linalg.lu_factor(A)\n",
    "x = linalg.lu_solve(LU, B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  6.  15. -23.]\n"
     ]
    }
   ],
   "source": [
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P=\n",
      " [[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n",
      "L=\n",
      " [[ 1.   0.   0. ]\n",
      " [ 0.5  1.   0. ]\n",
      " [ 0.5 -0.2  1. ]]\n",
      "U=\n",
      " [[ 2.   1.   1. ]\n",
      " [ 0.   2.5  1.5]\n",
      " [ 0.   0.  -0.2]]\n"
     ]
    }
   ],
   "source": [
    "import scipy\n",
    "\n",
    "P, L, U = scipy.linalg.lu(A)\n",
    "\n",
    "print('P=\\n', P)\n",
    "print('L=\\n', L)\n",
    "print('U=\\n', U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Cholesky decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "Cholesky decomposition with NumPy \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "\n",
    "A = np.array([\n",
    "    [10., -1., 2., 0.],\n",
    "    [-1., 11., -1., 3.],\n",
    "    [2., -1., 10., -1.],\n",
    "    [0., 3., -1., 8.]])\n",
    "B = np.array([6., 25., -11., 15.])\n",
    "\n",
    "L = np.linalg.cholesky(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.16227766  0.          0.          0.        ]\n",
      " [-0.31622777  3.3015148   0.          0.        ]\n",
      " [ 0.63245553 -0.24231301  3.08889696  0.        ]\n",
      " [ 0.          0.9086738  -0.25245792  2.6665665 ]]\n"
     ]
    }
   ],
   "source": [
    "print(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[10. -1.  2.  0.]\n",
      " [-1. 11. -1.  3.]\n",
      " [ 2. -1. 10. -1.]\n",
      " [ 0.  3. -1.  8.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.dot(L, L.T.conj()))  # A=L.L*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linalg.solve(L, B)  # L.L*.x=B; When L*.x=y, then L.y=B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linalg.solve(L.T.conj(), y)  # x=L*'.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.  2. -1.  1.]\n"
     ]
    }
   ],
   "source": [
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  6.]\n",
      " [ 25.]\n",
      " [-11.]\n",
      " [ 15.]]\n"
     ]
    }
   ],
   "source": [
    "print(np.mat(A) * np.mat(x).T)  # B=Ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The QR decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "QR decomposition with scipy \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import scipy.linalg as linalg\n",
    "\n",
    "\n",
    "A = np.array([\n",
    "    [2., 1., 1.],\n",
    "    [1., 3., 2.],\n",
    "    [1., 0., 0]])\n",
    "B = np.array([4., 5., 6.])\n",
    "\n",
    "Q, R = scipy.linalg.qr(A)  # QR decomposition\n",
    "y = np.dot(Q.T, B)  # Let y=Q'.B\n",
    "x = scipy.linalg.solve(R, y)  # Solve Rx=y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  6.  15. -23.]\n"
     ]
    }
   ],
   "source": [
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving with other matrix algebra methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Jacobi method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Solve Ax=B with the Jacobi method \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "\n",
    "def jacobi(A, B, n, tol=1e-10):\n",
    "    # Initializes x with zeroes with same shape and type as B\n",
    "    x = np.zeros_like(B)\n",
    "\n",
    "    for iter_count in range(n):\n",
    "        x_new = np.zeros_like(x)\n",
    "        for i in range(A.shape[0]):\n",
    "            s1 = np.dot(A[i, :i], x[:i])\n",
    "            s2 = np.dot(A[i, i + 1:], x[i + 1:])\n",
    "            x_new[i] = (B[i] - s1 - s2) / A[i, i]\n",
    "\n",
    "        if np.allclose(x, x_new, tol):\n",
    "            break\n",
    "\n",
    "        x = x_new\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([\n",
    "    [10., -1., 2., 0.], \n",
    "    [-1., 11., -1., 3.], \n",
    "    [2., -1., 10., -1.], \n",
    "    [0.0, 3., -1., 8.]])\n",
    "B = np.array([6., 25., -11., 15.])\n",
    "n = 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = [ 1.  2. -1.  1.]\n"
     ]
    }
   ],
   "source": [
    "x = jacobi(A, B, n)\n",
    "print('x', '=', x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Gauss-Seidel method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "Solve Ax=B with the Gauss-Seidel method \n",
    "\"\"\"\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def gauss(A, B, n, tol=1e-10):\n",
    "    L = np.tril(A)  # returns the lower triangular matrix of A\n",
    "    U = A-L  # decompose A = L + U\n",
    "    L_inv = np.linalg.inv(L)\n",
    "    x = np.zeros_like(B)\n",
    "\n",
    "    for i in range(n):\n",
    "        Ux = np.dot(U, x)\n",
    "        x_new = np.dot(L_inv, B - Ux)\n",
    "\n",
    "        if np.allclose(x, x_new, tol):\n",
    "            break\n",
    "\n",
    "        x = x_new\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([\n",
    "    [10., -1., 2., 0.], \n",
    "    [-1., 11., -1., 3.], \n",
    "    [2., -1., 10., -1.], \n",
    "    [0.0, 3., -1., 8.]])\n",
    "B = np.array([6., 25., -11., 15.])\n",
    "n = 100\n",
    "x = gauss(A, B, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = [ 1.  2. -1.  1.]\n"
     ]
    }
   ],
   "source": [
    "print('x', '=', x)"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}