Math629_midterm_coding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "142c7ea7",
   "metadata": {
    "id": "142c7ea7"
   },
   "source": [
    "**You are free to consult any source online, for syntax or conceptual understanding or for any other help you need but please write your own code**\n",
    "\n",
    "In Problem 2, don't worry if your results are bad or so, we are mainly looking at your approach to handle a new problem. Just put in your best effort and you will be fine, if you are using google colab, convert that into notebook and then submit it"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11ab29b2",
   "metadata": {
    "id": "11ab29b2"
   },
   "source": [
    "# Problem 1( 30 Points )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b22e2f02",
   "metadata": {
    "id": "b22e2f02"
   },
   "source": [
    "Till now in the course, we have looked at Pytorch for doing deep learning. However, it is important to know how neural networks are implemented from scratch because it helps a lot in doing debugging and it increases your fundamental understanding"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "400857c3",
   "metadata": {
    "id": "400857c3"
   },
   "source": [
    "In this problem, you will implement a very simple three layer neural network with Input $x$, $W_1, W_2, W_3$ be the weights and $b_1, b_2, b_3$ be the biases and $y$ be the output."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c420c8bd",
   "metadata": {
    "id": "c420c8bd"
   },
   "source": [
    "\n",
    "![Image](https://drive.google.com/uc?id=1TW85VcYzDNwTTWDKEOR1fSEAIw3Rtw5V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fc1764a",
   "metadata": {
    "id": "1fc1764a"
   },
   "source": [
    "We will work with BodyFat Data. Use 150 examples for training and the remainder  for  testing  with  mean  squared  error.  Activation  functions are rectified linear unit (ReLU) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94cc168a",
   "metadata": {
    "id": "94cc168a"
   },
   "source": [
    "You have to do the following:-\n",
    "\n",
    "a) Implement the forward pass <br>\n",
    "b) Implement the backward pass <br>\n",
    "c) Submit code and report on mean square error on training and testing data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88d66c3c",
   "metadata": {
    "id": "88d66c3c"
   },
   "source": [
    "# Notes:-"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5778a398",
   "metadata": {
    "id": "5778a398"
   },
   "source": [
    "1) seed : np.random.seed(0) <br>\n",
    "2) Initialize weights following a standard normal distribution <br>\n",
    "3) Initialize biases with zeros <br>\n",
    "3) First layer 64 hidden units, second layer 16 hidden units <br>\n",
    "4) **Use vectorized code** <br>\n",
    "5) Use gradient descent for training and use reasonable stopping criteria for termination, such as stop when loss doesn't decrease much, use your judgement <br>\n",
    "6) Learning rate: 1e-7 <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20b22df9",
   "metadata": {
    "id": "20b22df9"
   },
   "source": [
    "You can read the mat file using scipy functionality:-"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e56a60c-9902-4c57-9c8b-fc54a188a9c3",
   "metadata": {},
   "source": [
    "# Part 1 #\n",
    "\n",
    "*The discussion of the results is at the bottom of he code*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "id": "e5e47f9b",
   "metadata": {
    "id": "e5e47f9b"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/bv/j678dvtd34bc5mtymgsd2c6w0000gn/T/ipykernel_43857/2156593710.py:12: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  X_train = torch.tensor(torch.from_numpy(X_train).float(),dtype=torch.float64)\n",
      "/var/folders/bv/j678dvtd34bc5mtymgsd2c6w0000gn/T/ipykernel_43857/2156593710.py:13: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  X_test = torch.tensor(torch.from_numpy(X_test).float(),dtype=torch.float64)\n",
      "/var/folders/bv/j678dvtd34bc5mtymgsd2c6w0000gn/T/ipykernel_43857/2156593710.py:14: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  Y_train = torch.tensor(torch.from_numpy(Y_train).float(),dtype=torch.float64)\n",
      "/var/folders/bv/j678dvtd34bc5mtymgsd2c6w0000gn/T/ipykernel_43857/2156593710.py:15: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  Y_test = torch.tensor(torch.from_numpy(Y_test).float(),dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "import scipy.io as sio\n",
    "from sklearn.model_selection import train_test_split\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import numpy as np\n",
    "\n",
    "bodyfat = sio.loadmat('bodyfat_data.mat')\n",
    "X = bodyfat['X']\n",
    "Y = bodyfat['y']\n",
    "\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X,Y,test_size=98,shuffle=True)\n",
    "X_train = torch.tensor(torch.from_numpy(X_train).float(),dtype=torch.float64)\n",
    "X_test = torch.tensor(torch.from_numpy(X_test).float(),dtype=torch.float64)\n",
    "Y_train = torch.tensor(torch.from_numpy(Y_train).float(),dtype=torch.float64)\n",
    "Y_test = torch.tensor(torch.from_numpy(Y_test).float(),dtype=torch.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "id": "20c58d68-363a-45e5-9049-5839244eb27c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ScratchNeuralNet(nn.Module):\n",
    "    def __init__(self,inputSize, hiddenOne, hiddenTwo, outputSize):\n",
    "        super().__init__()\n",
    "        self.inputSize = inputSize\n",
    "        self.hiddenOne = hiddenOne\n",
    "        self.hiddenTwo = hiddenTwo\n",
    "        self.outputSize = outputSize\n",
    "        \n",
    "        torch.manual_seed(0)\n",
    "        \n",
    "        self.W1 = torch.normal(mean=torch.zeros(self.inputSize,self.hiddenOne),std=1/150*torch.ones(self.inputSize,self.hiddenOne)).double().requires_grad_()\n",
    "        self.b1 = torch.zeros(1,self.hiddenOne).double().requires_grad_()\n",
    "        \n",
    "        self.W2 = torch.normal(mean=torch.zeros(self.hiddenOne,self.hiddenTwo),std=1/150*torch.ones(self.hiddenOne, self.hiddenTwo)).double().requires_grad_()\n",
    "        self.b2 = torch.zeros(1,self.hiddenTwo).double().requires_grad_()\n",
    "        \n",
    "        self.W3 = torch.normal(mean=torch.zeros(self.hiddenTwo,self.outputSize),std=1/150*torch.ones(self.hiddenTwo,self.outputSize)).double().requires_grad_()\n",
    "        self.b3 = torch.zeros(1,self.outputSize).double().requires_grad_()\n",
    "        \n",
    "        self.relu = nn.ReLU()\n",
    "        \n",
    "    def forward(self,X):\n",
    "        a1 = torch.mm(X,self.W1) + self.b1\n",
    "        z1 = self.relu(a1)\n",
    "        \n",
    "        a2 = torch.mm(z1,self.W2) + self.b2\n",
    "        z2 = self.relu(a2)\n",
    "        \n",
    "        yhat = torch.mm(z2,self.W3) + self.b3\n",
    "        return yhat\n",
    "    \n",
    "    def cost(self, yhat, y):\n",
    "        AveMSE = torch.mean( (y-yhat)**2 ) \n",
    "        return AveMSE\n",
    "   \n",
    "    def zero_grads(self):\n",
    "        with torch.no_grad():\n",
    "            self.W1.grad.zero_()\n",
    "            self.W2.grad.zero_()\n",
    "            self.W3.grad.zero_()\n",
    "            self.b1.grad.zero_()\n",
    "            self.b2.grad.zero_()\n",
    "            self.b3.grad.zero_()\n",
    "\n",
    "    def updateParams(self, lr):\n",
    "        with torch.no_grad():\n",
    "            self.W1 -= lr*self.W1.grad\n",
    "            self.W2 -= lr*self.W2.grad\n",
    "            self.W3 -= lr*self.W3.grad\n",
    "            self.b1 -= lr*self.b1.grad\n",
    "            self.b2 -= lr*self.b2.grad\n",
    "            self.b3 -= lr*self.b3.grad\n",
    "        \n",
    "    \n",
    "    def backward(self,yhat,y):\n",
    "        loss = self.cost(yhat,y)\n",
    "        loss.backward()\n",
    "        self.updateParams(0.0001)\n",
    "        self.zero_grads()\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 404,
   "id": "e1b604ca-6cfa-4fb4-aab8-ac1f1f3b4a29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch:   0 Loss: tensor(432.4668, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   1 Loss: tensor(432.3133, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   2 Loss: tensor(432.1598, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   3 Loss: tensor(432.0061, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   4 Loss: tensor(431.8524, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   5 Loss: tensor(431.6986, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   6 Loss: tensor(431.5446, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   7 Loss: tensor(431.3905, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   8 Loss: tensor(431.2363, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   9 Loss: tensor(431.0818, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   10 Loss: tensor(430.9271, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   11 Loss: tensor(430.7722, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   12 Loss: tensor(430.6170, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   13 Loss: tensor(430.4615, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   14 Loss: tensor(430.3057, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   15 Loss: tensor(430.1495, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   16 Loss: tensor(429.9930, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   17 Loss: tensor(429.8359, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   18 Loss: tensor(429.6784, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   19 Loss: tensor(429.5203, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   20 Loss: tensor(429.3616, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   21 Loss: tensor(429.2023, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   22 Loss: tensor(429.0423, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   23 Loss: tensor(428.8815, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   24 Loss: tensor(428.7198, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   25 Loss: tensor(428.5571, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   26 Loss: tensor(428.3933, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   27 Loss: tensor(428.2284, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   28 Loss: tensor(428.0622, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   29 Loss: tensor(427.8945, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   30 Loss: tensor(427.7251, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   31 Loss: tensor(427.5540, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   32 Loss: tensor(427.3809, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   33 Loss: tensor(427.2056, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   34 Loss: tensor(427.0277, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   35 Loss: tensor(426.8463, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   36 Loss: tensor(426.6593, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   37 Loss: tensor(426.4686, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   38 Loss: tensor(426.2737, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   39 Loss: tensor(426.0740, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   40 Loss: tensor(425.8689, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   41 Loss: tensor(425.6575, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   42 Loss: tensor(425.4390, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   43 Loss: tensor(425.2123, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   44 Loss: tensor(424.9762, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   45 Loss: tensor(424.7290, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   46 Loss: tensor(424.4691, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   47 Loss: tensor(424.1944, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   48 Loss: tensor(423.9022, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   49 Loss: tensor(423.5896, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   50 Loss: tensor(423.2529, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   51 Loss: tensor(422.8873, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   52 Loss: tensor(422.4871, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   53 Loss: tensor(422.0455, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   54 Loss: tensor(421.5534, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   55 Loss: tensor(420.9997, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   56 Loss: tensor(420.3700, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   57 Loss: tensor(419.6462, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   58 Loss: tensor(418.8044, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   59 Loss: tensor(417.8136, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   60 Loss: tensor(416.6323, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   61 Loss: tensor(415.2049, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   62 Loss: tensor(413.4559, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   63 Loss: tensor(411.2811, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   64 Loss: tensor(408.5345, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   65 Loss: tensor(405.0091, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   66 Loss: tensor(400.4065, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   67 Loss: tensor(394.2906, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   68 Loss: tensor(386.0163, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   69 Loss: tensor(374.6226, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   70 Loss: tensor(358.6821, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   71 Loss: tensor(336.1264, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   72 Loss: tensor(304.1714, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   73 Loss: tensor(259.7819, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   74 Loss: tensor(201.7756, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   75 Loss: tensor(135.9070, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   76 Loss: tensor(79.5182, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   77 Loss: tensor(50.4128, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   78 Loss: tensor(43.5175, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "Epoch:   79 Loss: tensor(42.8808, dtype=torch.float64, grad_fn=<MeanBackward0>)\n",
      "\n",
      "Test:   79 Loss: tensor(59.0964, dtype=torch.float64, grad_fn=<MeanBackward0>)\n"
     ]
    }
   ],
   "source": [
    "model = ScratchNeuralNet(2,64,16,1)\n",
    "epochs = 80\n",
    "\n",
    "for epoch in range(epochs):\n",
    "    yhat = model(X_train)\n",
    "    model.backward(yhat,Y_train)\n",
    "    print(\"Epoch:\",\" \", epoch, \"Loss:\", model.cost(yhat,Y_train))\n",
    "\n",
    "\n",
    "yhatTest = model(X_test)\n",
    "model.backward(yhatTest, Y_test)\n",
    "print(\"\\nTest:\",\" \", epoch, \"Loss:\", model.cost(yhatTest,Y_test))\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "025cd91b-886e-49ba-9e03-468cc0ae0903",
   "metadata": {},
   "source": [
    "## Discussion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73f6c83e-d1b5-4fb1-a8fd-a6a4b07944c3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a62e493-430c-4246-9be4-26f238cd6900",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd151c82-cea0-4f9a-b237-91ea48b75dd9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c47883a8-54ae-4817-924f-206ef146c47e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "85313fca",
   "metadata": {
    "id": "85313fca"
   },
   "source": [
    "# Problem 2( 70 Points )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5502d743",
   "metadata": {
    "id": "5502d743"
   },
   "source": [
    "__Part 1 (Monthly Factor Data)__: Download the 5- factor Fama french factors monthly factors dataset from the below link:-\n",
    "    \n",
    "https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html\n",
    "\n",
    "Once you download the data, you will find that it contains a mix of annual and monthly factors and some texts at the top, so pd.read_csv will error out( try taking a look at csv file in text editor ). To get the monthly data, you can do the following:-\n",
    "\n",
    "df = pd.read_csv('F-F_Research_Data_5_Factors_2x3.csv', skiprows = 3).iloc[:-59,:]\n",
    "\n",
    "Use data from 1970-2010 as training data, 2010-2015 as validation data, 2015-2020 as a testing data. \n",
    "\n",
    "As you have seen from other courses, Factor models are heavily used in asset pricing as well as predicting future stock returns. For future stock return prediction, we typically have a seperate model for factor loadings. By factor loadings, we mean that coefficents of factor in the model. At time $t$ and by taking the dot product of 1) best estimation of factor loading for next time interval($t+1$) and 2) best estimate of factors at next time interval ($t+1$), one can come up with an estimation of future single stock return and create a long/short portfolio of stocks (long the stock with positive return prediction and short the ones that have negative predicted returns). However, most of these predictive models simply use the current value of the factors as the best guess for the factor value as time $t+1$. An alternative choice is using moving average of factor values as an estimate for the value of the factor at next time step.\n",
    "\n",
    "In this question, you are asked to analyze whether using a time series model (RNN/GRU/LSTM) would be able to give us better estimate for the value of the future factors values. To be more concrete, you need to try some recurrent neural networks and  use past values of the factors as input and the value of the factors at next time interval as the target. You can try different design for the network and its inputs, for example:\n",
    "\n",
    "- The input is 1-dimensional (the future value of each factor depends only on its own past).\n",
    "- The imput is 5-dimensional (the future value of each factor depends on history of all factors).\n",
    "\n",
    "You can either train 5 different neural networks to predict each of the 5 factors, or use one neural network to predict 5 factors at once\n",
    "\n",
    "Use data from 1970-2010 as training data, 2010-2015 as validation data, 2015-2020 as a testing data. You should tune the parameters using validation set only. After you have tuned the hyperparameters, use the model to make predictions on testing data and report the accuracy and performance of model. Compare the performance of the model with the two benchmarks mentioned above (last known value or moving average of previous values). Choose the window size for moving average as you desire\n",
    "\n",
    "__Part 2__ As you might have learnt from past courses, all these 5 factors are tradable portfolios by construction (for example SMB factor is a tradable portfolio loging small Caps and shorting Big Caps). Using the predictive model of previous step, design a trading strategy (for example, long the factor with highest estimated return and short the factor with smallest estimated return, or any other strategy of your choice). Calculate total return, sharpe ratio for the out of sample period.\n",
    "\n",
    "To be more precise, you can choose a starting cash amount such as 100,000 dollars, and see the performance of your trading strategy on this amount\n",
    "\n",
    "Suppose the factor predicted at time step is $F_{t+1}$, you can calculate the percentage return as $\\frac{F_{t+1} - F_{t}}{F_{t}}$. You can similarly calculate the return of true factors in testing data\n",
    "\n",
    "To be more precise, if the return vector of factors after prediction from neural network is $[0.1, 0.05, 0.07, 0.08, -0.09]$, then we will long the first factor and short the fifth factor. If the return vector of factors our testing dataset is $[0.02, 0.08, 0.54, 0.09, 0.05]$, then our PNL will be  $100,000*(-0.02 + 0.05)$ which is equal to $3000$\n",
    "\n",
    "You can sum this PNL across various time steps by calculating and that will give you the performace of trading strategy on entire testing data( out of sample data)\n",
    "\n",
    "As you can see, since the factors are tradables itself, we don't need to worry about factor loadings in calculating the PnL. You can experiment with any trading strategy of your choice based on these factors\n",
    "\n",
    "\n",
    "__Part 3__  Repeat the previous two steps by woking on daily factor data (you can download it from the same url) and perform a comparison between the two approaches\n",
    "\n",
    "When use download daily data, you will also find that it contains text in the beginning. Just set the parameter skiprows to 3 in pd.read_csv, and it will work fine and then you can proceed along similar lines\n",
    "\n",
    "Feel free to use your judgement on this problem and make assumptions anywhere you get stuck, just mention that in the notebook. The design of the problem is made open-ended intentionally\n",
    "\n",
    "__Feel free to use google colab for Keras or Pytorch.__<br>\n",
    "__Don't forget to use a lot of regularization to minimize the risk of overfitting__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1109,
   "id": "a09e068d-1653-4d84-b6e6-ecde007870a7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv('F-F_Research_Data_5_Factors_2x3.csv', skiprows = 3).iloc[:-59,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1110,
   "id": "b47f1315-2629-43fd-8791-b10cf7136745",
   "metadata": {},
   "outputs": [],
   "source": [
    "index1 = df[df[\"Unnamed: 0\"]==\"197001\"].index[0]\n",
    "index2 = df[df[\"Unnamed: 0\"]==\"200912\"].index[0]\n",
    "index3 = df[df[\"Unnamed: 0\"]==\"201512\"].index[0]\n",
    "index4 = df[df[\"Unnamed: 0\"]==\"202012\"].index[0]\n",
    "\n",
    "trainFeatures = df[[\"SMB\",\"HML\",\"RMW\",\"CMA\",\"RF\"]].loc[index1:index2]\n",
    "validationFeatures = df[[\"SMB\",\"HML\",\"RMW\",\"CMA\",\"RF\"]].loc[index2+1:index3]\n",
    "testFeatures = df[[\"SMB\",\"HML\",\"RMW\",\"CMA\",\"RF\"]].loc[index3+1:index4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1111,
   "id": "fd3cdb35-e263-45d8-9665-49eec400cd3c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train.shape =  (477, 2, 5)\n",
      "y_train.shape =  (477, 5)\n"
     ]
    }
   ],
   "source": [
    "def load_data(df, look_back):\n",
    "    data_raw = df.values # convert to numpy array\n",
    "    data = []\n",
    "    \n",
    "    # create all possible sequences of length look_back\n",
    "    for index in range(data_raw.shape[0] - look_back): \n",
    "        data.append(data_raw[index: index + look_back])\n",
    "    \n",
    "    data = np.array(data)\n",
    "    \n",
    "    x_train = data[:,:-1,:] #first (lookback - 1) elements\n",
    "    y_train = data[:,-1,:] #last element which is to be predicted based on x_train\n",
    "\n",
    "    return x_train, y_train\n",
    "\n",
    "look_back = 3 # sequence length\n",
    "x_train, y_train = load_data(trainFeatures, look_back)\n",
    "x_val, y_val = load_data(validationFeatures, look_back)\n",
    "x_test, y_test = load_data(testFeatures, look_back)\n",
    "print('x_train.shape = ',x_train.shape)\n",
    "print('y_train.shape = ',y_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1112,
   "id": "89342dda-07e4-4036-b27d-039e94e3fe1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = torch.from_numpy(np.asfarray(x_train)).type(torch.Tensor)\n",
    "x_test = torch.from_numpy(np.asfarray(x_test)).type(torch.Tensor)\n",
    "x_val = torch.from_numpy(np.asfarray(x_val)).type(torch.Tensor)\n",
    "y_train = torch.from_numpy(np.asfarray(y_train)).type(torch.Tensor)\n",
    "y_test = torch.from_numpy(np.asfarray(y_test)).type(torch.Tensor)\n",
    "y_val = torch.from_numpy(np.asfarray(y_val)).type(torch.Tensor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1113,
   "id": "885497f6-b358-4247-bc81-56ffc1a49bac",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dim = 5 # there is only single value at each node\n",
    "hidden_dim = 2\n",
    "num_layers = 2 #one LSTM stacked on top of each other\n",
    "output_dim = 5 #this is a regression problem now\n",
    "\n",
    "# Here we define our model as a class\n",
    "class LSTM(nn.Module):\n",
    "    def __init__(self, input_dim, hidden_dim, num_layers, output_dim):\n",
    "        super(LSTM, self).__init__()\n",
    "        # Hidden dimensions\n",
    "        self.hidden_dim = hidden_dim\n",
    "\n",
    "        # Number of hidden layers\n",
    "        self.num_layers = num_layers\n",
    "\n",
    "        self.n_outputs = output_dim\n",
    "\n",
    "        # batch_first=True causes input/output tensors to be of shape\n",
    "        # (batch_dim, seq_dim, feature_dim)\n",
    "        self.lstm = nn.LSTM(input_dim, hidden_dim, num_layers, batch_first=True,dropout = 0.5)\n",
    "\n",
    "        # Readout layer\n",
    "        self.fc = nn.Linear(hidden_dim, output_dim)\n",
    "\n",
    "    def forward(self, x):\n",
    "\n",
    "        self.batch_size = x.size(0)\n",
    "        # Initialize hidden state with zeros\n",
    "        h0 = torch.zeros(self.num_layers, self.batch_size, self.hidden_dim).requires_grad_()\n",
    "\n",
    "        # Initialize cell state\n",
    "        c0 = torch.zeros(self.num_layers, self.batch_size, self.hidden_dim).requires_grad_()\n",
    "\n",
    "        # We need to detach as we are doing truncated backpropagation through time (BPTT)\n",
    "        # If we don't, we'll backprop all the way to the start even after going through another batch\n",
    "\n",
    "        out, _ = self.lstm(x, (h0.detach(), c0.detach()))\n",
    "\n",
    "        final_out = self.fc(out[:,-1,:]) \n",
    "        final_out = final_out.view(-1,self.n_outputs)\n",
    "\n",
    "        return final_out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1114,
   "id": "4319c52d-c101-44a7-a425-78b7d700f531",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LSTM(\n",
      "  (lstm): LSTM(5, 2, num_layers=2, batch_first=True, dropout=0.5)\n",
      "  (fc): Linear(in_features=2, out_features=5, bias=True)\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "model = LSTM(input_dim=input_dim, hidden_dim=hidden_dim, output_dim=output_dim, num_layers=num_layers)\n",
    "\n",
    "loss = nn.MSELoss()\n",
    "\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.001)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1115,
   "id": "7f83e283-4e32-4b10-ab01-7d753abbc194",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Epoch  100 MSE:  5.94448709487915\n",
      "Epoch  100 Val MSE:  2.670475721359253\n",
      "\n",
      "Epoch  200 MSE:  5.874831199645996\n",
      "Epoch  200 Val MSE:  2.6509244441986084\n",
      "\n",
      "Epoch  300 MSE:  5.828480243682861\n",
      "Epoch  300 Val MSE:  2.66176700592041\n",
      "\n",
      "Epoch  400 MSE:  5.80929708480835\n",
      "Epoch  400 Val MSE:  2.6757209300994873\n"
     ]
    }
   ],
   "source": [
    "num_epochs = 500\n",
    "hist = np.zeros(num_epochs)\n",
    "hist2 = np.zeros(num_epochs)\n",
    "# Number of steps to unroll\n",
    "seq_dim =look_back-1  \n",
    "\n",
    "for t in range(num_epochs):\n",
    "    \n",
    "    # Forward pass\n",
    "    y_train_pred = model(x_train)\n",
    "\n",
    "    loss_val = loss(y_train_pred, y_train)\n",
    "\n",
    "    if t % 100 == 0 and t !=0:\n",
    "        print(\"\\nEpoch \", t, \"MSE: \", loss_val.item())\n",
    "        print(\"Epoch \",t,\"Val MSE: \",loss(model(x_val),y_val).item())\n",
    "\n",
    "    hist[t] = loss_val.item()\n",
    "    hist2[t] = loss(model(x_val),y_val).item()\n",
    "    # Zero out gradient, else they will accumulate between epochs\n",
    "    optimizer.zero_grad()\n",
    "\n",
    "    # Backward pass\n",
    "    loss_val.backward()\n",
    "\n",
    "    # Update parameters\n",
    "    optimizer.step()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1116,
   "id": "b2fa2c93-0649-44ad-b3b1-3cac5c3dc6c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0.5, 0, 'Number of Epochs'),\n",
       " Text(0, 0.5, 'Average Mean-Square Error'),\n",
       " Text(0.5, 1.0, 'Error vs. Num Epochs')]"
      ]
     },
     "execution_count": 1116,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA580lEQVR4nO3deXxU5dn4/881k8m+QCDsSwABNyRABDcU1FprEYviQm2V2qcuj3Whi33ax6326/Nr1Vprba3aVm3rWq1UrbbuikvVsCkoKEuAsCWE7OskuX5/3CdhiFkmITMhmev9es1rznqf68xyrnPuc859RFUxxhgTu3y9HYAxxpjeZYnAGGNinCUCY4yJcZYIjDEmxlkiMMaYGGeJwBhjYpwlAmNMm0TkIRH5f70dh4k8SwSmR4lIvojUiEhlyOue3o7rQHkbRRWRmSHDDhGRqNyIIyI3i0iw1edaGo1lm/7PEoGJhDNVNTXk9d22JhKRuDaG+buyoK5Of4D2Ar25h/xEq891QC/GYvoRSwQmakRksYi8IyK/EpFi4GZvT/teEXlBRKqAuSJymIi8ISKlIrJWROaHlPGF6Vst43wRyWs1bImIPOt1nyEin4hIhYhsF5EfdGEVHgaOEpGT2lm/fBE5NaT/ZhH5q9ed7R1RfEtEtolIiYhcLiJHi8hH3rp2+8jJK/tqEdkkIntE5HYR8XnjfCJyvYhsEZFCEfmziGSEzHuCiLzrxbBNRBaHFD1QRP7pfV7vi8gEbx7xvsdCESkXkY9F5Mjuxm96lyUCE22zgE3AUOBWb9jXve404H3gOeAlYAhwFfCIiEwOKSN0+rdblf8cMFlEJraa/lGv+4/AZaqaBhwJvNaF2KuB/wuJuztmAROB84G7gP8FTgWOAM5rL8mEaQGQC0wHzgIu8YYv9l5zgfFAKnAPgIiMBV4EfgNkATnAqpAyLwB+CgwENrBv3U8DTgQmARnAeUDxAcRuepElAhMJS729y+bXd0LG7VDV36hqg6rWeMP+oarvqGoTbkOUCvxcVetV9TXgeWBRSBkt06tqbeiCVbUa+Efz9F5COBR41pskCBwuIumqWqKqK7q4bvcBY0TkK12cr9nPVLVWVV8CqoDHVLVQVbcDy4BpHcx7XqvP9fVW43+hqntVdSsuyTR/ZhcCd6rqJlWtBH4MXOBVzX0deEVVH1PVoKoWq+qqkDKfUdUPVLUBeAT3/YD7HNNwn62o6qequrN7H4npbZYITCR8TVUHhLweCBm3rY3pQ4eNALZ5SaHZFmBkJ2WEepR9G8GvA0u9BAFwDnAGsEVE3hSRYztbmVCqWgf8zHt1x+6Q7po2+lM7mPfJVp/r3FbjQz+XLbjPEu99S6txcbijstHAxg6WuSuku7o5Pi9B3wP8FigUkftFJL2DcsxBzBKBiba2rrIJHbYDGN1cv+0ZA2zvpIxQLwNZIpKDSwjN1UKo6oeqehau2mkp8GTYke/zIDAAOLvV8CogOaR/WDfKPhCjQ7rH4D5LvPexrcY14JLQNmBCdxamqner6gzgcFwV0Q+7U47pfZYIzMHmfdye53UiEhCROcCZwOPhFqCqQeBvwO1AJi4xICLxInKhiGR405QDTe2X1G75DcBNwI9ajVqFq3IJiEgusLCrZR+gH4rIQBEZDVwDPOENfwxYIiLjRCQVd57jiZDqnlNF5DwRiRORQV4C7ZB3knuWiARwCbCWbnyW5uBgicBEwnOy//Xuz4Q7o6rW4zb8XwH2AL8DLlLVdV2M4VHcSdi/eRu8Zt8E8kWkHLgcV3+OiIzxYh0TZvmPAa3rxG/A7V2X4E6wPtp6pgN0fqvPtVJEhoSM/wewHJeQ/ok7MQ7wJ+AvwFvAZtxG+yoA73zCGcD3cZfHrgKmhhFLOvAAbl234E4U334gK2d6j9iDaYzp+8Td2DZRVTf0diym77EjAmOMiXGWCIwxJsZZ1ZAxxsQ4OyIwxpgY94VGvw52gwcP1uzs7N4Owxhj+pTly5fvUdWstsb1uUSQnZ1NXl5e5xMaY4xpISJb2htnVUPGGBPjLBEYY0yMs0RgjDExzhKBMcbEOEsExhgT4ywRGGNMjLNEYIwxMS5mEkFRRR3/7/lPWL2tlDXby6isa+h8JmOMiQERvaFMRAYAf8A9JFyBS1T1vZDxAvwa1x56NbC4G8+QDct7m4p58N18/vD2ZgDi43yMHphEQ5Ny2LB0hmUkMmFIKuMHp1BV18DMcZkkxftJiPNHIhxjjDloRPrO4l8D/1LVhSISz/6P8QP38JGJ3msWcK/33uPmTx3BUSMzWLerAlVl5bZSVm8rJTUhjpXbSti9tu4L88TH+Zg8NI09lXWMG5xCVloCx08YzMShqYzOTGZQSjwulxljTN8VsdZHRSQD97Sj8drOQkTkPuANVX3M618PzFHV1k9+apGbm6s93cREY5PS0NTE9pIa8vJLSE+KY3VBGfUNTazdUUZGUoA128uprm+gpDrYMl9KvJ+0xAADkgMADMtIZHBqAtmDkhk7KIVBqfEEG5WxmcmMHZSMiNDU5D4Kn88SiDEmekRkuarmtjUukkcE44Ai4EERmYp7hN41qloVMs1I3MOzmxV4w/ZLBCJyKXApwJgx4T5JMHx+n+D3+Rmflcr4rFQATj9y+BemU1U2FFaypbiabSXVbCmupqwmyN6qekqq69lQWMnaHeU8VfHFo4vBqfEMTU9kd3kd9Q2N5IwZyPQxA0iI83PosDTSk+Iorqxn5MAkRmcmk54Y6PH1NMaYtkQyEcQB04GrVPV9Efk18D+457p2iareD9wP7oigR6PsAhFh4tA0Jg5N63C62mAjGworKa8NEufzsaGwkhVbSyipqicp4Ke0JsinO8t567OiNueP9/tIivczZWQGmSnxiIBfhOzBKWQkBUiK91Nd18Cs8YOorm9kQlYKqQlxxPlj5ty/MaYHRTIRFAAFqvq+1/8ULhGE2g6MDukf5Q3r0xIDfo4cmdHSP3NcJl+ftf+RTFOTUlhRR1K8n892V7CrrJbBqQms3VHGxqIqgo1NrN5WyraSauqCTYjA31e2/9EkBnwMS08kLTGAovhEOGJEBqMzkxiUEs+A5Hhqg40MS09keEYSPp87EhqalmjVVMbEuIglAlXdJSLbRGSyqq4HTgE+aTXZs8B3ReRx3Eniso7OD/QnPp8wLCMRgKOzM1uGHzthULvzFFbU4hOhqq4BVXh7w56W/m0l1eytqmdbSQ1pCQGaVHnh452U1QTbLQ8gIynA4NR40pMCpCcGvPc4kgJ+BqbEMyA5wKCUBBLifKQnxZGaECAp4GfMoNbn/Y0xfVWkrxq6CnjEu2JoE/AtEbkcQFV/D7yAu3R0A+7y0W9FOJ4+bUiaSxyDUxMAyB6c0uk8zSe4S6rqCfh97C6vZXd5LapQ29DIpzsrKK8JUl4bpLS6nq17qymvCVJV30BtsKnNMn0CwzOSAEhNiGNQajxNqgxKSWBQqjvymDgkjcyUeOL8wqiBSWSmJJCaEEdiwEdVXSMDkgMkBuzSXGMOBn3umcWRuGrItG1PZR3BxiaKK+upa2iivCZIUUUd63dXsKW4isSAn2BjE7vL66isa6C0OkhlXbDdBNJaVpq7wgpgaHoiowYms6eyjtSEOGrqG5k6egAAU0dnMHJAEvUNTaQkxJEc77fLdo3pot66asj0cc1HHs17/+EINjahCjvLamhSqKlvZGdZDeW1QSpr3VFGcoKf4sp6CkqqWb+rgvg4Hx8VlPH8RzvJSApQVddAUsDPE3nb2lyGT1zi8IlQVd9AemKAI0akMyQtgWEZSYi4S4JFIDHOT1K8nxEDksgZPYCqugaKK+upqA0ybcxA4uN8+ARLLCamWSIwPSrgXbk0dtC+aqvDR6SHNW9dQyPxIVc+bSyqoqqugY1FlRRV1JGcEEdVXQNVdQ1sKa4mzif4fcKu8lre+qyIqvrGLsXq9wmNTcqYzGSy0hKoqmvA7xNGD0xmWEYijU3KsIxEMpICfL67gqNGDSA1MY7BqfEMTk0gMeBnS3E1k4amogpNqgxMjreT76bPsURgDhqtm/M4ZIi7p6O5iqgjNfWNiLgjkoKSGnwijB2UTF2wiepgA+9sKGZvVR0Bv49BqQmoKq+tK2RAUoDNxdU0NjWRGIinpr6RFVtLqK5vJOCX/W4ghHYf+drC7xOyByVTWh1kxtiBBOJ8JAX8BPw+4v3CzrJaJg1NY8SAJAoragn4fZw0KYs9lXVsLKoiMeDjkKxUaoKN7ubEwSmkxPvZtKeK8YNTeuTIpbFJaVJtSdrG2DkCYzpQUFJNUUUdhw5LZ/3uCuJ8QlFlHXsq6iirCTIkPZGVW0vISksg4POxs6yWrXurSYr383FBKY2q7CqrJTHOj+LuPt9YVEm4fzsRSE8MUFYTJCngp1GV7EHJJAb87oglPdFdhhzwMzQjkcamJrJSE8hMSaCiNsikYWnUBRsZn5XK7vJaCivq+FveNoZnJPGXb8+kJtjI5j1VxPt9BBuVitogIuKdh4HJQ9PYW13PpqIqZo3L9D6TGrLSEmhSpaFJ27z5salJUVxiVFVU7W763tbROQJLBMZEWVl1kJpgI3F+4aF38hGBEQOSGD84hZEDk/i8sJLymiB1wSZ2ldeyflcFCQEf6YkB4uN8LN9SQlVdA8MzEtlSXM2IAUmU1tRTUhVEVSnxyu/MwOQApTXBTpNSUsBPTbCRgckBqusbqWtoIiMpQE2wEVVl5rhM4nw+PtlZTlVdA2MykymprifO5+PkQ4ewYmsJtcFG5k4ewqDUBMYOSuaDzXv515pdJCf4Oe3wYQxOjWdIuktkpdVBBiQHmD0xiysfWcFhw9PJHpTMtDEDaVJ3j8yO0hre37yXJV+axJ7KOjKT49lYVMnKraXMmzqctMQAH+bvZcLgVDKSAqQlxuHzCX9YtomMpABnTx+F30tMZdVBNu6pJCMpQEKcj6y0BBLi/KgqNcFGkuPbrjipqW8kMeCjSd3VeWmJAYoq6oj3+0hNjGspv3kb29vnoSwRGBNDmpqUbSXVJMT5efnT3aQnxpGZEs+ApHjSk+II+H08/F4+O0prmTgklQlZqTR5Gz0UAnFCwd4afvnyZwCcdvhQMlPieeXT3ZTXNPClI4ZSXhNk1MBk6hoa+TB/L5W1DcwaN4j0pDh2lNayuqCUSu9+l7TEOCpqv9jse/M5mvaI0GmS8gk06b739kwdlcGWvdWUelV9WWnucuaiirovNEk/cUgqYzKT+SB/L1V1DSycMYoNhZXsKK0le3AyJxwymPc2FfPB5r2cOXUEO0tr+Xh7GWdMGcaTeQUAzBqXyYJpI3kibxsrt5YyaWgqV508kQ/z97KjtIaxg1IYPTCJtMQA1fUNjM9K5bgJg1izvZxNeyr5ZGc5Ewan8uUjh/HU8gLKqus55bChHDUqo9sJxRKBMabLSqvryUgK7LfhUdUubYiq6hqIj/ORv6eKZZ/vYcG0kawqKGXl1lKuOvkQBKiobaAm2EhBSQ2NTcr4rBTe/KyIlVtLOHPqCI4cmUFpVZDfvr6BwWnxJAX8TBk1gNpgI8+t3sH0MQMpqa4nKd7PyAFJXPP4KgB+fUEOu8treWdDMcu3lFBZ18CQtAR+csZhPPbBVtKTAuwsc+eTymuC5BdXM3JAEgG/0KjK8RMGs720huVb3DkjcEdRJdVBMlPi2VtV35Mfd4eak+K3TxjHDfMO72YZlgiMMTHinQ17KCip5vyj9zXr0tik+ATqGpravJGxNthIwO9rqc4Jpao0KQhug1xYUcfg1AR8Ap/sLCfg95Ec72drcTXDMhIZMSCJFVtLeH1dIefmjmZ7SQ3DByTy6qeFXHxcNhsLK7n3jY1cdNxYauobGZqeyO3/Xs+bnxUxPiuFkyZl8ej7W5mQlUpDUxN3npfD6MxkXvx4J5OGpTF9zMBufS6WCIwx5iC3tbiazNR4UhPiCDY2EfD7unwE1hG7ocwYYw5yoe13NV/aG60TzHYhsTHGxDhLBMYYE+MsERhjTIyzRGCMMTHOEoExxsQ4SwTGGBPjLBEYY0yMs0RgjDExzhKBMcbEOEsExhgT4yLaxISI5AMVQCPQ0LqdCxGZA/wD2OwN+ruq3hLJmIwxxuwvGm0NzVXVPR2MX6aq86IQhzHGmDZY1ZAxxsS4DhOBiPhFZMkBlK/ASyKyXEQubWeaY0VktYi8KCJHtBPHpSKSJyJ5RUVFBxCOMcaY1jpMBKraCCw6gPJPUNXpwFeAK0XkxFbjVwBjVXUq8BtgaTtx3K+quaqam5WVdQDhGGOMaS2cqqF3ROQeEZktItObX+EUrqrbvfdC4BlgZqvx5apa6XW/AAREZHDXVsEYY8yBCOdkcY73Hno1jwIndzSTiKQAPlWt8LpPa1UGIjIM2K2qKiIzcYmpOMzYjTHG9IBOE4Gqzu1m2UOBZ7wn7MQBj6rqv0Tkcq/c3wMLgStEpAGoAS7QvvbsTGOM6eM6TQQikgHcBDTX778J3KKqZR3Np6qbgKltDP99SPc9wD1dCdgYY0zPCuccwZ9wN4Wd573KgQcjGZQxxpjoCeccwQRVPSek/6cisipC8RhjjImycI4IakTkhOYeETkeV59vjDGmHwjniOBy4M/euQKAEuDiyIVkjDEmmjpMBCLiB76pqlNFJB3ctf9RicwYY0xUdJgIVLWxuVrIEoAxxvRP4VQNrRSRZ4G/AVXNA1X17xGLyhhjTNSEkwgScXf7ht5JrIAlAmOM6QfCOUdQrKo/iFI8xhhjoiyc1kePj1IsxhhjekE4VUOr7ByBMcb0X3aOwBhjYlw4rY9+KxqBGGOM6R3tniMQkSdDun/RatxLkQzKGGNM9HR0snhiSPeXWo2z50UaY0w/0VEi6OgBMfbwGGOM6Sc6OkeQLCLTcMkiyesW75UUjeCMMcZEXkeJYCdwp9e9K6S7ud8YY0w/0G4iOIBnFRtjjOlDwnkwjTHGmH4soolARPJF5GMRWSUieW2MFxG5W0Q2iMhHIjI9kvEYY4z5onDuLD5Qc1V1TzvjvoK7THUiMAu413s3xhgTJZ0eEXh77d8QkRu9/jEiMrOHln8W8Gd1/gMMEJHhPVS2McaYMIRTNfQ74FhgkddfAfw2zPIVeElElovIpW2MHwlsC+kv8IbtR0QuFZE8EckrKioKc9HGGGPCEU4imKWqVwK1AKpaAsSHWf4JqjodVwV0pYic2J0gVfV+Vc1V1dysLLup2RhjelI4iSDoPaBGAUQkC2gKp3BV3e69FwLPAK2rlLYDo0P6R3nDjDHGREk4ieBu3EZ8iIjcCrwN/F9nM4lIioikNXcDpwFrWk32LHCRdx7iGKBMVXd2ZQWMMcYcmM4eVekDNgPXAafgmpf4mqp+GkbZQ4FnRKR5OY+q6r9E5HIAVf098AJwBrABqAasyWtjjImyDhOBqjaJyG9VdRqwrisFq+omYGobw38f0q3AlV0p1xhjTM8Kp2roVRE5R7xde2OMMf1LOIngMtzziutEpFxEKkSkPMJxGWOMiZJwHlWZFo1AjDHG9I6wmpgQkYG4ZiASm4ep6luRCsoYY0z0dJoIROS/gGtw1/ivAo4B3gNOjmhkxhhjoiKccwTXAEcDW7xnFEwDSiMZlDHGmOgJJxHUqmotgIgkqOo6YHJkwzLGGBMt4ZwjKBCRAcBS4GURKQG2RDIoY4wx0RPOVUMLvM6bReR1IAP4V0SjMsYYEzXhnCweE9K72XsfBmyNSETGGGOiKpyqoX/iWh4V3OWj44D1wBERjMsYY0yUhFM1NCW033uu8H9HLCJjjDFR1eWH16vqCuy5wsYY02+Ec47geyG9PmA6sCNiERljjImqcM4RhLY11IA7Z/B0ZMIxxhgTbeGcI/hpNAIxxhjTO8KpGnoO73nFbVHV+T0akTGmVwWDQQoKCqitre3tUEw3JCYmMmrUKAKBQNjzhFM1tAl338Bfvf5FwG7cncbGmH6moKCAtLQ0srOzsedR9S2qSnFxMQUFBYwbNy7s+cJJBMeram5I/3MikqeqS7ocpTHmoFdbW2tJoI8SEQYNGkRRUVGX5gvn8tEUERkfsqBxQEoX4zPG9CGWBPqu7nx34SSCJcAbIvKGiLwJvI5rmjrcoPwislJEnm9j3GIRKRKRVd7rv8IP3RjTHxUXF5OTk0NOTg7Dhg1j5MiRLf319fUdzpuXl8fVV1/d6TKOO+64Hon1jTfeICMjoyW+nJwcXnnllR4pO5rCuWroXyIyETjUG7ROVeu6sIxrgE+B9HbGP6Gq3+1CecaYfmzQoEGsWrUKgJtvvpnU1FR+8IMftIxvaGggLq7tTVdubi65ubltjgv17rvv9kisALNnz+b557+wn9tCVVFVfD5fm/3t6Wg9e1q7kYjI0SIyDMDb8E8FbgFuF5HMcAoXkVHAV4E/9ECsxpgYtXjxYi6//HJmzZrFddddxwcffMCxxx7LtGnTOO6441i/fj3g9tDnzZsHuCRyySWXMGfOHMaPH8/dd9/dUl5qamrL9HPmzGHhwoUceuihXHjhhai6iyRfeOEFDj30UGbMmMHVV1/dUm448vPzmTx5MhdddBFHHnkky5Yt269/27Zt/PCHP+TII49kypQpPPHEEy3xzJ49m/nz53P44Yf3yGcXjo7SzX3AqQAiciLwc+AqIAe4H1gYRvl3Adex/01prZ3jlf8ZsERVt7WeQEQuBS4FGDNmTOvRxpgI+elza/lkR3mPlnn4iHRuOrPrbVYWFBTw7rvv4vf7KS8vZ9myZcTFxfHKK6/wk5/8hKef/uJ9ruvWreP111+noqKCyZMnc8UVV3zhssqVK1eydu1aRowYwfHHH88777xDbm4ul112GW+99Rbjxo1j0aJF7ca1bNkycnJyWvqffvpp/H4/n3/+OQ8//DDHHHMM+fn5+/U//fTTrFq1itWrV7Nnzx6OPvpoTjzxRABWrFjBmjVrunTVz4HqKBH4VXWv130+cL+qPg08LSKrOitYROYBhaq6XETmtDPZc8BjqlonIpcBD9PGs5BV9X5c8iE3N7fdexqMMf3Xueeei9/vB6CsrIyLL76Yzz//HBEhGAy2Oc9Xv/pVEhISSEhIYMiQIezevZtRo0btN83MmTNbhuXk5JCfn09qairjx49v2RgvWrSI+++/v81ltFU1lJ+fz9ixYznmmGNahoX2v/322yxatAi/38/QoUM56aST+PDDD0lPT2fmzJlRTQLQSSIQkThVbQBOwdsjD2O+ZscD80XkDFzz1eki8ldV/UbzBKpaHDL9H4Dbwg/dGBNp3dlzj5SUlH0XK95www3MnTuXZ555hvz8fObMmdPmPAkJCS3dfr+fhoaGbk1zoPG21R/ufNHQ0dmKx4A3ReQfQA2wDEBEDgHKOitYVX+sqqNUNRu4AHgtNAl4ZQ0P6Z2PO6lsjDEdKisrY+TIkQA89NBDPV7+5MmT2bRpE/n5+QAtdfg9Zfbs2TzxxBM0NjZSVFTEW2+9xcyZM3t0GV3RbiJQ1VuB7wMPASdo8xkUN89V3V2giNwiIs3NUlwtImtFZDVwNbC4u+UaY2LHddddx49//GOmTZvWY3vwoZKSkvjd737H6aefzowZM0hLSyMjI6PNaZvPETS/nnrqqU7LX7BgAUcddRRTp07l5JNP5rbbbmPYsGE9vRphk33b9zAmFrnUq6/vNbm5uZqXl9ebIRjTr3366accdthhvR1Gr6usrCQ1NRVV5corr2TixIksWdI3GlRo6zsUkeWtWolo0dUH01ze3cCMMaYveeCBB8jJyeGII46grKyMyy67rLdDipiu3q1g950bY2LCkiVL+swRwIHq6hHBmRGJwhhjTK8J53kECcA5QDYQ19ygkareEtHIjDHGREU4VUP/wF0uuhzoShtDxhhj+oBwEsEoVT094pEYY4zpFeGcI3hXRKZEPBJjjAHmzp3Lv//97/2G3XXXXVxxxRXtzjNnzhyaLys/44wzKC0t/cI0N998M3fccUeHy166dCmffPJJS/+NN97YI81KH+zNVYdzRHACsFhENuOqhgRQVT0qopEZY2LSokWLePzxx/nyl7/cMuzxxx/nttvCa4HmhRde6Payly5dyrx581pa/rzllp47FXowN1cdzhHBV4CJwGm4q4bmYVcPGWMiZOHChfzzn/9seQhNfn4+O3bsYPbs2VxxxRXk5uZyxBFHcNNNN7U5f3Z2Nnv27AHg1ltvZdKkSZxwwgktTVWDu0fg6KOPZurUqZxzzjlUV1fz7rvv8uyzz/LDH/6QnJwcNm7cyOLFi1vuFH711VeZNm0aU6ZM4ZJLLqGurq5leTfddBPTp09nypQprFu3Lux1PViaqw7nwTRbAERkCK7xOGNMrLj2WvAeEtNjcnLgrrvaHZ2ZmcnMmTN58cUXOeuss3j88cc577zzEBFuvfVWMjMzaWxs5JRTTuGjjz7iqKParpxYvnw5jz/+OKtWraKhoYHp06czY8YMAM4++2y+853vAHD99dfzxz/+kauuuor58+czb948Fi7cv5X92tpaFi9ezKuvvsqkSZO46KKLuPfee7n22msBGDx4MCtWrOB3v/sdd9xxB3/4wxcfwXIwN1fd6RGBiMwXkc+BzcCbQD7w4gEv2Rhj2tFcPQSuWqj5eQBPPvkk06dPZ9q0aaxdu3a/+vzWli1bxoIFC0hOTiY9PZ358+e3jFuzZg2zZ89mypQpPPLII6xdu7bDeNavX8+4ceOYNGkSABdffDFvvfVWy/izzz4bgBkzZrQ0VNfa7NmzWbVqVctrwoQJAN1qrhro0eaqw6lY+hlwDPCKqk4TkbnANzqZxxjTH3Sw5x5JZ511FkuWLGHFihVUV1czY8YMNm/ezB133MGHH37IwIEDWbx4MbW1td0qf/HixSxdupSpU6fy0EMP8cYbbxxQvM1NWXenGeuDobnqcM4RBL3nBvhExKeqrwOdPxTUGGO6KTU1lblz53LJJZe0HA2Ul5eTkpJCRkYGu3fv5sUXO66YOPHEE1m6dCk1NTVUVFTw3HPPtYyrqKhg+PDhBINBHnnkkZbhaWlpVFRUfKGsyZMnk5+fz4YNGwD4y1/+wkknndQTq9qhaDVXHc4RQamIpOKeR/CIiBQCVT0eiTHGhFi0aBELFixoqSKaOnUq06ZN49BDD2X06NEcf/zxHc4/ffp0zj//fKZOncqQIUM4+uijW8b97Gc/Y9asWWRlZTFr1qyWjf8FF1zAd77zHe6+++79mpNOTEzkwQcf5Nxzz6WhoYGjjz6ayy/vWhucrc8RXH/99eTmdrxPvWDBAt577z2mTp2KiLQ0V92VE9Lh6LQZahFJwT2YxgdcCGQAj7R6uljUWDPUxkSWNUPd93W1GepwrhqqEpGxwERVfVhEkgF/j0RrjDGm14Vz1dB3gKeA+7xBI4GlEYzJGGNMFIVzsvhK3IPoywFU9XNgSCSDMsYYEz3hJII6Va1v7hGROCD851saY/qcrjzC1hxcuvPdhZMI3hSRnwBJIvIl4G/Ac53MY4zpoxITEykuLrZk0AepKsXFxSQmdq0RiHAuH/0f4NvAx8BlwAvAF++fboeI+IE8YLuqzms1LgH4MzADKAbOV9X8cMs2xvS8UaNGUVBQQFFRUW+HYrohMTGRUaNGdWmecK4aagIe8F7dcQ3wKZDexrhvAyWqeoiIXAD8Aji/m8sxxvSAQCDQY00XmL6h3UQgIh91NGM4zVCLyCjgq8CtwPfamOQs4Gav+yngHhERtWNSY4yJmo6OCJpwJ4UfxZ0TqOlG+XcB1wFp7YwfCWwDUNUGESkDBgF7QicSkUuBSwHGjBnTjTCMMca0p92TxaqaAywCUnHJ4FbgCFxd/5bOChaReUChqi4/0CBV9X5VzVXV3KysrAMtzhhjTIgOrxpS1XWqepOqTscdFfwZWBJm2ccD80UkH3gcOFlE/tpqmu3AaGi5LDUDd9LYGGNMlHSYCERkpIh8X0TexjU9vQS4N5yCVfXHqjpKVbOBC4DXVLV189XPAhd73Qu9aez8gDHGRFFHJ4vfxNXtPwl8i3176vEikqmqe7uzQBG5BchT1WeBPwJ/EZENwF5cwjDGGBNF7bY+6lXpNI8Mnaj54fXjIxta26z1UWOM6bputT7qVekYY4zp58JpYsIYY0w/ZonAGGNinCUCY4yJcWElAhE5QUS+5XVniYg1RGKMMf1EOE8ouwn4EfBjb1AAaH1jmDHGmD4qnCOCBcB8oApAVXfQfttBxhhj+phwEkG9d7evAohISmRDMsYYE03hJIInReQ+YID3IPtX6P6zCYwxxhxkwnkwzR3eIyrLgcnAjar6csQjM8YYExXhPKoSb8NvG39jjOmHOk0EIlLB/m0NAZThnkP8fVXdFInAjDHGREc4RwR3AQW4h9MIroXQCcAK4E/AnAjFZowxJgrCOVk8X1XvU9UKVS1X1fuBL6vqE8DACMdnjDEmwsJJBNUicp6I+LzXeUCtN84eImOMMX1cOIngQuCbQCGw2+v+hogkAd+NYGzGGGOiIJzLRzcBZ7Yz+u2eDccYY0y0hXPVUCLwbeAIILF5uKpeEsG4jDHGREk4VUN/AYYBXwbeBEYBFZEMyhhjTPSEkwgOUdUbgCpVfRj4KjArsmEZY4yJlnASQdB7LxWRI4EMYEhnM4lIooh8ICKrRWStiPy0jWkWi0iRiKzyXv/VtfCNMcYcqHBuKLtfRAYC1wPPAqnADWHMVwecrKqVIhIA3haRF1X1P62me0JV7eojY4zpJR0mAhHxAeWqWgK8BYwPt2Cv6epKrzfgvey+A2OMOch0WDWkqk3Add0tXET8IrIKdw/Cy6r6fhuTnSMiH4nIUyIyup1yLhWRPBHJKyoq6m44xhhj2hDOOYJXROQHIjJaRDKbX+EUrqqNqpqDu9JopneOIdRzQLaqHoVr3fThdsq5X1VzVTU3KysrnEUbY4wJUzjnCM733q8MGaZ0rZqoVEReB04H1oQMLw6Z7A/AbeGWaYwxpmeEc2fxuO4ULCJZQNBLAknAl4BftJpmuKru9HrnA592Z1nGGGO6L5w7i5OB7wFjVPVSEZkITFbV5zuZdTjwsIj4cVVQT6rq8yJyC5Cnqs8CV4vIfKAB2AssPoB1McYY0w3iLu7pYAKRJ4DlwEWqeqSXGN716v6jLjc3V/Py8npj0cYY02eJyHJVzW1rXDgniyeo6m14N5apajXuATXGGGP6gXASQb1Xx68AIjIBd7OYMcaYfiCcq4ZuBv4FjBaRR4Djsbp8Y4zpN8K5auglEVkOHIOrErpGVfdEPDJjjDFREc5VQ8/hHlz/rKpWRT4kY4wx0RTOOYI7gNnAJ14zEAu9h9UYY4zpB8KpGnoTeNO7H+Bk4DvAn4D0CMdmjDEmCsI5WYx31dCZuOYmptNOm0DGGGP6nnDOETwJzMRdOXQP8KbXKqkxxph+IJwjgj8Ci1S1EUBEThCRRap6ZSfzGWOM6QPCOUfwbxGZJiKLgPOAzcDfIx6ZMcaYqGg3EYjIJGCR99oDPIFrm2hulGIzxhgTBR0dEawDlgHzVHUDgIgsiUpUxhhjoqaj+wjOBnYCr4vIAyJyCtbYnDHG9DvtJgJVXaqqFwCHAq8D1wJDROReETktSvEZY4yJsE7vLFbVKlV9VFXPxD17eCXwo4hHZowxJirCaWKihaqWeA+SPyVSARljjImuLiUCY4wx/Y8lAmOMiXGWCIwxJsZFLBGISKKIfCAiq0VkrYj8tI1pEkTkCRHZICLvi0h2pOIxxhjTtkgeEdQBJ6vqVCAHOF1Ejmk1zbeBElU9BPgV8IsIxmOMMaYNEUsE6lR6vQHvpa0mO4t9TVo/BZwiInbTmjHGRFFEzxGIiF9EVgGFwMuq+n6rSUYC2wBUtQEoAwa1Uc6lIpInInlFRUWRDNkYY2JORBOBqjaqag7uRrSZInJkN8u5X1VzVTU3KyurR2M0xphYF5WrhlS1FNdMxemtRm0HRgOISByQARRHIyZjjDFOJK8ayhKRAV53EvAlXIumoZ4FLva6FwKvqWrr8wjGGGMiKKxnFnfTcOBh76H3PuBJVX1eRG4B8lT1WdzTz/4iIhuAvcAFEYzHGGNMGyKWCFT1I2BaG8NvDOmuBc6NVAzGGGM6Z3cWG2NMjLNEYIwxMc4SgTHGxDhLBMYYE+MsERhjTIyzRGCMMTHOEoExxsQ4SwTGGBPjLBEYY0yMi2QTEweX3bvhgQegsBAOPxwaG2HlSggG4eKLob4e/H6YMAF27IChQyEjA4YM6e3IjTHd0djo/tM9LT8f0tJg0BdazO9cUxP4Dr7979hJBG+8ATfcAAkJUFfnhqWnQ3k5/PnP7c83cqRLCIceCmvXui8xLQ2OPx4GD3b9mZnwyScwdy7s2QMFBTBlihv3la+4ZQaD7oeZmBiV1TXmoPPpp+6/MnSo6//wQ1i/Hr7xDdff0OD+M21tKBsa4Ktfhf/+bzjrLDds1y5QheHDvzj9Qw/BkiWwaRMMHLhveE2N+w+GPv/q44/hwQfh9ts7TxwNDTBunNuZXLsWysrguuvcDuZVV8HJJ+8/fVkZ3HknfP/7sHEjTJ8Ob70Fs2e3XX5ZmdvW3HknLFrkhlVXQ1KS26ZcfDH88Icdx9gN0tca+8zNzdW8vLyuzxgMwubNbo9/1y43bMQI2LrVfUFJSVBVBR99BGPHwt69UFTkfqyVlbBtm8vmw4e7I4YNG8Jf9uDBroxg0O1FJCW5ZVRUQG4ufPCBOyIJBOCoo9wP7Oij4ZRT3DRjxrjk89Zb7n3KFJg2zSW0mhqXqMrKYNiwyOwBmf6peQNzIA8F3LvXbViTkzuerqnJ/TYHDIBf/ALmzXM7WeD+gykpbiNZUwMXXghf+5rbiJ96Krz0Eowf73a0AO66C95+G/79b/f/uPFG+M1voKQEHn0UTjwR5sxx/9E77nD/sZoat3F9+WX43vfcNC+95LYBt9/u/j+HHQa1tXDrre797rth4kT4+c+htNRt/P/1L1iwoP31XLAA4uLgoovg3XdhyxYX06mnuviWL3fT/e538PDDcMYZLv6VK934U07ZV9Ytt7gd2Hfegfvug8WL3XxXXNGNLwpEZLmq5rY5LmYSQU8rL3d7+LW17s+QleWOLCZMcD/YTZtc4vngA9i5023A6+tdEqqrcz/g2lqXhGbMcHtKu3bB++/D6NEuCVVVdRyDiNsjahYX5/Z+0tIgPn7/hDNhglv+tm0wc6Zb3p49MGmSm8fncwll5043fVycm87nc+Vs2OD2VHbtcknI73d7R3FxB7YhMW1T/eLnunWr2xFpbIQVK1wVxbBhbqO2fj0cccT+0xcUwN//DqNGud9iYiJMngyrVrnv8Otfd2Wdfrr7PsvL3W+gosId4U6fDtde6zaQf/2r2yDNnr1vY/XSS3D11S4JjB3r4jjvPFizxm24duxwG/hDDoFXX21/XU87zZXVurs/CQTcjmBbUlI6/683++gjtyPYDZYI+pLmDUBNjTuUHjjQ/eGLimDWLLehXrnS7VmMHev+bIWFbo8jP9/tfRxyiDu3UVgI773n9qYKClwiOOww92NqajqwOAcOdHswiYnuUH/YMPdeV+diTU93e38lJW59Jk50e6CZmS62jAyXlLZvd0c8gcC+hLpundtwLV7shicluaSVm+uS54YN7nMqL4czz3RJd80a9wdJSnLlhP6uf/ITtwc6erQrb/du+NKX3Ofy9tvus/j4YxfH2LEupsmTITvbLa+uzi0/Odnt4X34oVufhgb4/HM37ccfu+9s7ly3Ed29221c165130NaGnz5y27d1q1zy/b73V5jba3buB5yCLz+uit37Vq3x3z44a6syZPhqadcgq+vdxvwZpMmwWefuWnnzHHLWr3aHUFWV3f+Xfp8X/w9nHCC+2xgX3Vq85FnqBEj4Jhj3Lq/+qqLrbWhQ906tPb1r7vP5pVXXP9ll8G997qqjz/9ycVVXAxTp7oYzjvPfX8/+5n7LpOT3e9l1Sq3d/3ss/DLX7qyRo+Gc85xRw/NBg6E116DX/3KxXT77W743Llw+eWuWueii+DFF13V74gRbrn19e6z3bTJfSbXXgtPP+3W+X//F447Dv75Txf/nXe6eC+80P1WwK3f8ce7o5MbbtgXT2Ki++5D3Xab+02PGOF+xytXuqOYe+91v5f6+m6fY7BEEMtC9yybu7dvdxuIwYPdH3vPHveD3LnTbUgnTXLTPfec22jX1rrpCwvdBmHbNvcHyc52w3ftcn/03bv3VZ+Vl7vpRdyGafdu98PfuXPfOZpmbW2I+prWe3VJSW5D0ZbZs93GSxXOPReeecZVPYQ65ZR9e9FjxrijgeZ5/X446SSXLDZudJ/x+PHuu3znnX0b/+xsePJJ99lv3uyOAnbtcjsD//mPS4Zjxrjv88EH3V7/7be772vuXLfBfewxSE113V/+sov7xRfdd7l9O1x6qUtS4L7vF190R59ZWXDsse4CjbPOghdecNUaTz/tvv8XXnCJYPly+NGP4P/+zyWUZjU17je0Y4eLO1TrE66hJ4U/+sj9rtPT3dHqL38J11/vPs9HHtl3fgJc0j/jDLeRbdZcn79yJeTktP39taf1Udzq1bB0qav6aR5eU+N+J3v2uM/1Bz9wSaKxcd9RV0LC/uW++66b5kc/ctVU3WSJwBxcSkrcYXJiovsDx8e7P31iotsgrl7tqjkaG93wqiq3N9ecmE491e0pFha6vbTx493G7fPPXbVGdrbbUDQ2ug3oiSe6o5Tmk4QZGW6DtHu3q+4aNsztfa9e7aYZNcrt9RcUuGkDAbfRqK93G6FTTnHL3rrVxbl2rduIlZS4DYjf7/p//3u3ZzhggNtI1Ne798zMfUkiKcntOSYkuGSyevW++Xfv3rfhWrfObbRTUjr+bINBlwhuuw2++U23frFu+3aXHFpvYNvS/JsZNy7ycXXF8uUuIR7AFUeWCIwxJsZ1lAgOvgtajTHGRJUlAmOMiXGWCIwxJsZFLBGIyGgReV1EPhGRtSJyTRvTzBGRMhFZ5b1ubKssY4wxkRPJJiYagO+r6goRSQOWi8jLqvpJq+mWqeq8CMZhjDGmAxE7IlDVnaq6wuuuAD4FRkZqecYYY7onKucIRCQbmAa838boY0VktYi8KCJHtDEeEblURPJEJK+oqCiSoRpjTMyJeCIQkVTgaeBaVS1vNXoFMFZVpwK/AZa2VYaq3q+quaqam5WVFdF4jTEm1kT0hjIRCQDPA/9W1TvDmD4fyFXVPR1MUwRs6WZIg4F2y+6nbJ1jg61zbDiQdR6rqm3uSUfsZLGICPBH4NP2koCIDAN2q6qKyEzcEUpxR+W2tyJhxpTX3p11/ZWtc2ywdY4NkVrnSF41dDzwTeBjEVnlDfsJMAZAVX8PLASuEJEGoAa4QPtamxfGGNPHRSwRqOrbQIcN1avqPcA9kYrBGGNM52LtzuL7ezuAXmDrHBtsnWNDRNa5z7U+aowxpmfF2hGBMcaYViwRGGNMjIuZRCAip4vIehHZICL/09vx9BQR+ZOIFIrImpBhmSLysoh87r0P9IaLiNztfQYficj03ou8+9pr0LA/r7eIJIrIB95d+GtF5Kfe8HEi8r63bk+ISLw3PMHr3+CNz+7VFegmEfGLyEoRed7r79frC+5+KhH52GuIM88bFtHfdkwkAhHxA78FvgIcDiwSkcN7N6oe8xBweqth/wO8qqoTgVe9fnDrP9F7XQrcS9/U3KDh4cAxwJXe99mf17sOONm7Cz8HOF1EjgF+AfxKVQ8BSoBve9N/Gyjxhv/Km64vugbXTlmz/r6+zeaqak7IPQOR/W2rar9/Acfi7m5u7v8x8OPejqsH1y8bWBPSvx4Y7nUPB9Z73fcBi9qari+/gH8AX4qV9QaScc2zzMLdZRrnDW/5nQP/Bo71uuO86aS3Y+/ieo7yNnon41ookP68viHrnQ8MbjUsor/tmDgiwLV6ui2kv4D+3RLqUFXd6XXvArwnoPe/z6FVg4b9er29apJVQCHwMrARKFXVBm+S0PVqWWdvfBkwKKoBH7i7gOuAJq9/EP17fZsp8JKILBeRS71hEf1tR/LOYnMQUFUVkX55jXDrBg1dqyZOf1xvVW0EckRkAPAMcGjvRhQ5IjIPKFTV5SIyp5fDibYTVHW7iAwBXhaRdaEjI/HbjpUjgu3A6JD+Ud6w/mq3iAwH8N4LveH95nPwGjR8GnhEVf/uDe736w2gqqXA67iqkQEi0rxDF7peLevsjc+gk3a8DjLHA/O9higfx1UP/Zr+u74tVHW7916IS/gzifBvO1YSwYfARO+Kg3jgAuDZXo4pkp4FLva6L8bVoTcPv8i70uAYoCzkcLPPEGm3QcN+u94ikuUdCSAiSbhzIp/iEsJCb7LW69z8WSwEXlOvErkvUNUfq+ooVc3G/V9fU9UL6afr20xEUsQ90RERSQFOA9YQ6d92b58YieIJmDOAz3D1qv/b2/H04Ho9BuwEgrj6wW/j6kZfBT4HXgEyvWkFd/XURuBjXJPfvb4O3VjnE3D1qB8Bq7zXGf15vYGjgJXeOq8BbvSGjwc+ADYAfwMSvOGJXv8Gb/z43l6HA1j3OcDzsbC+3vqt9l5rm7dVkf5tWxMTxhgT42KlasgYY0w7LBEYY0yMs0RgjDExzhKBMcbEOEsExhgT4ywRmIOWiKiI/DKk/wcicnMPlf2QiCzsfMoDXs65IvKpiLzeani2iNR4LUw2vy7qweXOaW6x05jOWBMT5mBWB5wtIv+fqu7p7WCaiUic7mvvpjPfBr6j7hnerW1U1Zyei8yY7rEjAnMwa8A9o3VJ6xGt9+hFpNJ7nyMib4rIP0Rkk4j8XEQuFNeW/8ciMiGkmFNFJE9EPvPatmlu2O12EfnQa9/9spByl4nIs8AnbcSzyCt/jYj8wht2I+7mtz+KyO3hrrSIVIrIr8Q9d+BVEcnyhueIyH+8uJ4JaZP+EBF5RdyzClaErGOqiDwlIutE5BHvjmy8z+QTr5w7wo3L9GO9fSedvezV3guoBNJxzfJmAD8AbvbGPQQsDJ3We58DlOKa6k3AtbvyU2/cNcBdIfP/C7czNBF3V3Yirk33671pEoA8YJxXbhUwro04RwBbgSzcUfZrwNe8cW/Qxt2euKbDa9h3Z/QqYLY3ToELve4bgXu87o+Ak7zuW0LW5X1ggdediGumeg6uBc5R3jq+h0tKg3BNFTffTDqgt79ne/X+y44IzEFNVcuBPwNXd2G2D1V1p6rW4W69f8kb/jFuA9zsSVVtUtXPgU241jxPw7Xdsgq3gR2ESxQAH6jq5jaWdzTwhqoWqasyegQ4MYw4N6p7+Ejza5k3vAl4wuv+K3CCiGTgNtpvesMfBk702qUZqarPAKhqrapWh8RboKpNuESTjUsOtbijlLOB5mlNDLNEYPqCu3B17Skhwxrwfr8i4gPiQ8bVhXQ3hfQ3sf95sdbtqyiu7ZarQjbO41S1OZFUHchKHIDutgMT+jk04h7o0oBrzfIpYB7uqMjEOEsE5qCnqnuBJ9n3WEJw1UUzvO75QKAbRZ8rIj6vTn08rsrk38AVXjPXiMgkrxXIjnwAnCQig8U9FnUR8GYn83TEx74WNr8OvK2qZUCJiMz2hn8TeFNVK4ACEfmaF2+CiCS3V7C4ZzhkqOoLuHMvUw8gTtNP2FVDpq/4JfDdkP4HgH+IyGrcXm139ta34jbi6cDlqlorIn/AVaGs8E6uFgFf66gQVd0pIv+DayJZgH+q6j86msczwauCavYnVb0bty4zReR6XLvz53vjLwZ+723oNwHf8oZ/E7hPRG7BtUJ7bgfLTMN9bolerN8LI07Tz1nro8YcZESkUlVTezsOEzusasgYY2KcHREYY0yMsyMCY4yJcZYIjDEmxlkiMMaYGGeJwBhjYpwlAmOMiXH/P/SRXtdZY79IAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "\n",
    "hist = pd.DataFrame(hist)#\n",
    "hist.rename(columns={0:\"Training Error\"},inplace=True)\n",
    "hist2 = pd.DataFrame(hist2)\n",
    "hist2.rename(columns={0:\"Validation Error\"},inplace = True)\n",
    "\n",
    "fig = sns.lineplot(data = hist)\n",
    "fig = sns.lineplot(data = hist2,palette=['r'])\n",
    "fig.set(xlabel=\"Number of Epochs\",ylabel=\"Average Mean-Square Error\",title = \"Error vs. Num Epochs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55943fb6-ed28-419c-a9d7-204c6946eed2",
   "metadata": {},
   "source": [
    "I now consider my model fully optimized. my methods for regulating over fitting are setting the dropout probability to 50%. I also have reduced the number of hidden dimension to just 2 and the sequence length to 1. I am genuinely suprised by the model parameters, but the larger my model tended to get in any of these parameter, the more quickly I saw evidence of severe overfitting. Now I will compare this LSTM model to baseline methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1117,
   "id": "fa4d18e4-d082-4ca4-9b93-0afa389d6b3a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The model has loss of 5.2334394454956055 over the TEST SET\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f\"\\nThe model has loss of {loss(model(x_test),y_test).item()} over the TEST SET\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1118,
   "id": "46a8b422-1749-401e-aae4-82f858b3f573",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The loss on the TEST SET from using the factors at time t as a proxy for the factors at time t+1 gives us an average MSE of 8.808276176452637\n",
      "\n"
     ]
    }
   ],
   "source": [
    "lastVal = x_test[:,-1,:]\n",
    "lastSeqPredictionLoss = loss(lastVal,y_test).item()\n",
    "print(f'\\nThe loss on the TEST SET from using the factors at time t as a proxy for the factors at time t+1 gives us an average MSE of {lastSeqPredictionLoss}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1128,
   "id": "cccb322e-2dae-426e-a6c7-e1d7c9366dc0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The loss on the TEST SET from using a moving average over the last 5 periods as a proxy for the factors at time t+1 gives us an average MSE of 6.325165668558522\n",
      "\n"
     ]
    }
   ],
   "source": [
    "temp_X_Test = np.array(x_test[:,0,:])\n",
    "MALength = 5\n",
    "MA = np.zeros((temp_X_Test.shape[0]-MALength,5))\n",
    "\n",
    "for i in range(temp_X_Test.shape[0]-MALength):\n",
    "    mySlice = temp_X_Test[i:i+MALength]\n",
    "    MA[i,:] = np.average(mySlice,axis=0)\n",
    "\n",
    "temp = torch.from_numpy(MA)\n",
    "MAPredLoss = loss(temp,y_test[MALength:]).item()\n",
    "print(f'\\nThe loss on the TEST SET from using a moving average over the last 5 periods as a proxy for the factors at time t+1 gives us an average MSE of {MAPredLoss}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dbd249d-4af1-4482-b729-14e4a0b6daaa",
   "metadata": {},
   "source": [
    "Overall, we can see that our simple LSTM network has improved the performance quite considerably from both baseline methods. Given the simplicity of the model I believe this absolutely warrants the use of this model over the baseline methods.\n",
    "\n",
    " \n",
    " \n",
    "\n",
    "I spent a considerable amount of time Tweaking the following to get my results:\n",
    "\n",
    "Learning Rate: Adjusted much lower to reduce the instability of Error\n",
    "\n",
    "Hidden Dimension: I realized that reducing the size of the hidden dimension greatly helped with reducing overfitting \n",
    "\n",
    "Num Epochs: I observed many different epoch sizes and found a suitable stopping time based off my architecture and learning rate\n",
    "\n",
    "Sequence Size: The sequence size was reduced to 2. This seems odd, but I did not find the results that I wanted with larger sequence sizes\n",
    "\n",
    "Dropout Rate: I used a dropout rate of the 50% which seems quite high, but lower values for dropout suggested that there was still considerable overfitting to address"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a7afff3-21d0-4d2f-8f23-accde00abeb6",
   "metadata": {},
   "source": [
    "# Part 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06315afb-d490-427b-a044-66d362767a48",
   "metadata": {},
   "source": [
    "I am going to try a long/short strategy for this section by following the recommended path above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1077,
   "id": "a14e4fc1-66fd-4741-886e-ea4d24200cb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictedFactors = model(x_test)\n",
    "predictedReturns = ((predictedFactors - x_test[:,-1,:])/x_test[:,-1,:]).detach().numpy()\n",
    "actualReturns = ((y_test - x_test[:,-1,:])/x_test[:,-1,:]).detach().numpy()\n",
    "actualReturns[np.ma.masked_invalid(actualReturns).mask] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1078,
   "id": "f401675e-341b-428a-82ca-df486927aeb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "strategy = np.zeros(predictedReturns.shape)\n",
    "strategy[np.max(predictedReturns,axis=1).reshape(strategy.shape[0],1)==predictedReturns] = 1\n",
    "strategy[np.min(predictedReturns,axis=1).reshape(strategy.shape[0],1)==predictedReturns] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1079,
   "id": "f70ad9e0-d9e0-4012-b761-4b0f3e4afc93",
   "metadata": {},
   "outputs": [],
   "source": [
    "Returns = actualReturns*strategy\n",
    "Returns = np.sum(Returns,axis=1).reshape(Returns.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1080,
   "id": "feadaa0d-8338-4704-a7fe-423b4cc8eb04",
   "metadata": {},
   "outputs": [],
   "source": [
    "CashAmount = 100000\n",
    "Returns[np.ma.masked_invalid(Returns).mask] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1163,
   "id": "886462f1-998a-4d8e-bf65-166e2687c5f4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Total PNL sums up to be: $-825525.3905430436\n"
     ]
    }
   ],
   "source": [
    "PNL = np.sum(Returns*CashAmount)\n",
    "print(f\"The Total PNL sums up to be: ${PNL}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8965bb83-f8c1-4aaf-ba78-73f3f40dd398",
   "metadata": {},
   "source": [
    "We can see that the returns are not a particularly great way to create a strategy. This is actually because calculated returns for values that move around zero produce highly unstable results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34142de6-3893-44d9-b4fe-13e4dd10fd75",
   "metadata": {},
   "source": [
    "As a second way to evaluate this, let us consider buying 10,000 shares for the factor that has the greatest difference from t to t+1 and we can short 10,000 shares that has the most negative difference from t to t+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1164,
   "id": "3f3d1867-98b3-493e-a98b-aff84edda37b",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictedFactors = model(x_test)\n",
    "predictedDiff = ((predictedFactors - x_test[:,-1,:])).detach().numpy()\n",
    "actualDiff = ((y_test - x_test[:,-1,:])).detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1165,
   "id": "b21856b6-d555-4c2a-9868-264b5a48306b",
   "metadata": {},
   "outputs": [],
   "source": [
    "strategy = np.zeros(predictedDiff.shape)\n",
    "strategy[np.max(predictedDiff,axis=1).reshape(strategy.shape[0],1)==predictedDiff] = 1\n",
    "strategy[np.min(predictedDiff,axis=1).reshape(strategy.shape[0],1)==predictedDiff] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1166,
   "id": "aeaad629-6145-41b6-88de-6b3d1cc5aa02",
   "metadata": {},
   "outputs": [],
   "source": [
    "Diff1 = actualDiff*strategy\n",
    "Diff = np.sum(Diff1,axis=1).reshape(Diff1.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1177,
   "id": "98669cd1-22b6-4d06-b854-a81ff5a3fb2c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The total gain over the test set for buying 10,000 shares of the factor with the biggest expected difference and shorting 10,000 shares with the smallest expected difference is $2150899.988543242\n"
     ]
    }
   ],
   "source": [
    "totalNumberShares = 10000\n",
    "PNL = Diff*totalNumberShares\n",
    "print(f'The total gain over the test set for buying 10,000 shares of the factor with the biggest expected difference and shorting 10,000 shares with the smallest expected difference is ${np.sum(PNL)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ce2629-1047-4bf6-ab4d-1514e05f200f",
   "metadata": {},
   "source": [
    "We can see that this is a much better strategy given the values that the Fama/French factors take on. It shows that we can reliably predict the direction of movement in the Fama/French factors and that we can implement a profitable strategy with this knowledge. Given that the maximum value of the Fama/French factors in the test set is about 8.2, we can expect that our strategy will never cost us more than $82,000 in one month"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e08bf2a4-57d0-49ba-a471-dc5f53958d22",
   "metadata": {},
   "source": [
    "Also note that because I have decided to look at the differences as apposed to the returns (For reasons stated above), I cannot calculated the sharpe ratio. This is primarily because my initial capital requirements are sometimes negative, thus suggesting that I require no capital for some month to month trades. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e77be6bb-9cd8-4922-ad6e-8c06c2c69993",
   "metadata": {},
   "source": [
    "# Part 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1208,
   "id": "ac8d398d-82b0-4afa-8e52-9fa8ee3c0a96",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "df = pd.read_csv('F-F_Research_Data_5_Factors_2x3_daily.csv', skiprows = 3).iloc[:-59,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1209,
   "id": "75b5d4cd-fde4-44bd-8dbe-7483141390fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "index1 = df[df[\"Unnamed: 0\"]==19700102].index[0]\n",
    "index2 = df[df[\"Unnamed: 0\"]==20091231].index[0]\n",
    "index3 = df[df[\"Unnamed: 0\"]==20151231].index[0]\n",
    "index4 = df[df[\"Unnamed: 0\"]==20201231].index[0]\n",
    "\n",
    "trainFeatures = df[[\"SMB\",\"HML\",\"RMW\",\"CMA\",\"RF\"]].loc[index1:index2]\n",
    "validationFeatures = df[[\"SMB\",\"HML\",\"RMW\",\"CMA\",\"RF\"]].loc[index2+1:index3]\n",
    "testFeatures = df[[\"SMB\",\"HML\",\"RMW\",\"CMA\",\"RF\"]].loc[index3+1:index4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1216,
   "id": "3db588c2-cceb-4aaf-9b6e-be165da49829",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_train.shape =  (10067, 29, 5)\n",
      "y_train.shape =  (10067, 5)\n"
     ]
    }
   ],
   "source": [
    "def load_data(df, look_back):\n",
    "    data_raw = df.values # convert to numpy array\n",
    "    data = []\n",
    "    \n",
    "    # create all possible sequences of length look_back\n",
    "    for index in range(data_raw.shape[0] - look_back): \n",
    "        data.append(data_raw[index: index + look_back])\n",
    "    \n",
    "    data = np.array(data)\n",
    "    \n",
    "    x_train = data[:,:-1,:] #first (lookback - 1) elements\n",
    "    y_train = data[:,-1,:] #last element which is to be predicted based on x_train\n",
    "\n",
    "    return x_train, y_train\n",
    "\n",
    "look_back = 30 # sequence length\n",
    "x_train, y_train = load_data(trainFeatures, look_back)\n",
    "x_val, y_val = load_data(validationFeatures, look_back)\n",
    "x_test, y_test = load_data(testFeatures, look_back)\n",
    "print('x_train.shape = ',x_train.shape)\n",
    "print('y_train.shape = ',y_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1217,
   "id": "cb5d1729-996a-4cc6-9b6f-c4d02ce3b72b",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = torch.from_numpy(np.asfarray(x_train)).type(torch.Tensor)\n",
    "x_test = torch.from_numpy(np.asfarray(x_test)).type(torch.Tensor)\n",
    "x_val = torch.from_numpy(np.asfarray(x_val)).type(torch.Tensor)\n",
    "y_train = torch.from_numpy(np.asfarray(y_train)).type(torch.Tensor)\n",
    "y_test = torch.from_numpy(np.asfarray(y_test)).type(torch.Tensor)\n",
    "y_val = torch.from_numpy(np.asfarray(y_val)).type(torch.Tensor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1259,
   "id": "a3d4fc29-9002-48bb-b753-4d6a6005f899",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_dim = 5 # there is only single value at each node\n",
    "hidden_dim = 5\n",
    "num_layers = 2 #one LSTM stacked on top of each other\n",
    "output_dim = 5 #this is a regression problem now\n",
    "\n",
    "# Here we define our model as a class\n",
    "class LSTM(nn.Module):\n",
    "    def __init__(self, input_dim, hidden_dim, num_layers, output_dim):\n",
    "        super(LSTM, self).__init__()\n",
    "        # Hidden dimensions\n",
    "        self.hidden_dim = hidden_dim\n",
    "\n",
    "        # Number of hidden layers\n",
    "        self.num_layers = num_layers\n",
    "\n",
    "        self.n_outputs = output_dim\n",
    "\n",
    "        # batch_first=True causes input/output tensors to be of shape\n",
    "        # (batch_dim, seq_dim, feature_dim)\n",
    "        self.lstm = nn.LSTM(input_dim, hidden_dim, num_layers, batch_first=True,dropout = 0.1)\n",
    "\n",
    "        # Readout layer\n",
    "        self.fc = nn.Linear(hidden_dim, output_dim)\n",
    "\n",
    "    def forward(self, x):\n",
    "\n",
    "        self.batch_size = x.size(0)\n",
    "        # Initialize hidden state with zeros\n",
    "        h0 = torch.zeros(self.num_layers, self.batch_size, self.hidden_dim).requires_grad_()\n",
    "\n",
    "        # Initialize cell state\n",
    "        c0 = torch.zeros(self.num_layers, self.batch_size, self.hidden_dim).requires_grad_()\n",
    "\n",
    "        # We need to detach as we are doing truncated backpropagation through time (BPTT)\n",
    "        # If we don't, we'll backprop all the way to the start even after going through another batch\n",
    "\n",
    "        out, _ = self.lstm(x, (h0.detach(), c0.detach()))\n",
    "\n",
    "        final_out = self.fc(out[:,-1,:]) \n",
    "        final_out = final_out.view(-1,self.n_outputs)\n",
    "\n",
    "        return final_out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1294,
   "id": "bf86cf52-3751-417e-ae3a-1e6cf3a9be9c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LSTM(\n",
      "  (lstm): LSTM(5, 5, num_layers=2, batch_first=True, dropout=0.1)\n",
      "  (fc): Linear(in_features=5, out_features=5, bias=True)\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "model = LSTM(input_dim=input_dim, hidden_dim=hidden_dim, output_dim=output_dim, num_layers=num_layers)\n",
    "\n",
    "loss = nn.MSELoss()\n",
    "\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.003)\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1295,
   "id": "e797e60b-8fa5-4669-a647-a7fe480a19ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Epoch  20 MSE:  0.19090870022773743\n",
      "Epoch  20 Val MSE:  0.1500094085931778\n",
      "\n",
      "Epoch  40 MSE:  0.17739947140216827\n",
      "Epoch  40 Val MSE:  0.13609562814235687\n"
     ]
    }
   ],
   "source": [
    "num_epochs = 55\n",
    "hist = np.zeros(num_epochs)\n",
    "hist2 = np.zeros(num_epochs)\n",
    "# Number of steps to unroll\n",
    "seq_dim =look_back-1  \n",
    "\n",
    "for t in range(num_epochs):\n",
    "    \n",
    "    # Forward pass\n",
    "    y_train_pred = model(x_train)\n",
    "\n",
    "    loss_val = loss(y_train_pred, y_train)\n",
    "\n",
    "    if t % 20 == 0 and t !=0:\n",
    "        print(\"\\nEpoch \", t, \"MSE: \", loss_val.item())\n",
    "        print(\"Epoch \",t,\"Val MSE: \",loss(model(x_val),y_val).item())\n",
    "\n",
    "    hist[t] = loss_val.item()\n",
    "    hist2[t] = loss(model(x_val),y_val).item()\n",
    "    # Zero out gradient, else they will accumulate between epochs\n",
    "    optimizer.zero_grad()\n",
    "\n",
    "    # Backward pass\n",
    "    loss_val.backward()\n",
    "\n",
    "    # Update parameters\n",
    "    optimizer.step()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1296,
   "id": "3c4dccf9-4d17-4140-a6df-a87789e42569",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0.5, 0, 'Number of Epochs'),\n",
       " Text(0, 0.5, 'Average Mean-Square Error'),\n",
       " Text(0.5, 1.0, 'Error vs. Num Epochs')]"
      ]
     },
     "execution_count": 1296,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABBT0lEQVR4nO3dd3hUZfbA8e9JSCUhISGA1NA7BAhNur2CBVGsrK7dVXHVdX/rqusu23R3XdeyYl9FsYKoWBEEBekBpEkLEFB6D4GU8/vjvYEhhmRIZjIp5/M895m59869c+4Q5sx9q6gqxhhjTFFhoQ7AGGNM5WQJwhhjTLEsQRhjjCmWJQhjjDHFsgRhjDGmWJYgjDHGFMsShDHmpIjIKyLyp1DHYYLPEoSpECKSKSKHROSAz/JUqOMqL+/LUkWkt8+21iJSIR2MROQREckt8rnuqYj3NtWfJQhTkS5U1Tif5Y7iXiQitYrZFn4yb3Syry+nXUAof1G/VeRzTQxhLKYasQRhQk5ERovItyLyLxHZCTzi/TJ/VkSmiMhBYKiIdBCR6SKyR0SWicgwn3P87PVF3uNyEZlfZNsYEZnsPT9PRJaLyH4R2Swi957EJbwKdBWRwSe4vkwROcNn/RERed17nurdgfxCRDaJyG4RuUVEeonIEu9ay3yn5Z37ThFZJyI7ROQxEQnz9oWJyIMiskFEtonI/0QkwefYASIyy4thk4iM9jl1XRH52Pu85ohIK+8Y8f4dt4nIPhFZKiKdyxq/CS1LEKay6AOsAxoAY71tV3rP44E5wIfA50B94FfAeBFp53MO39d/U+T8HwLtRKRNkde/4T1/EbhZVeOBzsBXJxF7NvBnn7jLog/QBrgceAL4HXAG0AkYeaLk46eLgXSgBzAcuN7bPtpbhgItgTjgKQARaQ58AvwHSAHSgAyfc14B/AGoC6zh2LWfBQwC2gIJwEhgZzliNyFkCcJUpEner9HC5UaffVtU9T+qmqeqh7xtH6jqt6pagPuCigP+qqpHVPUr4CNglM85jr5eVXN831hVs4EPCl/vJYr2wGTvJblARxGpo6q7VXXhSV7bc0AzETn3JI8r9EdVzVHVz4GDwJuquk1VNwMzge4lHDuyyOc6rcj+v6nqLlXdiEs+hZ/ZVcA/VXWdqh4Afgtc4RXxXQl8qapvqmququ5U1Qyfc05U1bmqmgeMx/37gPsc43GfrajqClX9sWwfiQk1SxCmIl2kqok+y/M++zYV83rfbY2ATV6yKLQBaFzKOXy9wbEvxyuBSV7iALgUOA/YICJfi0i/0i7Gl6oeBv7oLWWx1ef5oWLW40o49u0in+vQIvt9P5cNuM8S73FDkX21cHdxTYG1JbznTz7Pswvj8xL3U8DTwDYRGScidUo4j6nELEGYyqK4Vj++27YATQvLzz3NgM2lnMPXF0CKiKThEkVh8RKqOk9Vh+OKryYBb/sd+TEvA4nAJUW2HwRifdYbluHc5dHU53kz3GeJ99i8yL48XHLaBLQqy5up6pOq2hPoiCtquq8s5zGhZwnCVBVzcL9U7xeRCBEZAlwITPD3BKqaC7wDPAYk4RIGIhIpIleJSIL3mn1AwYnPdMLz5wEPA78psisDV3QTISLpwIiTPXc53ScidUWkKXAX8Ja3/U1gjIi0EJE4XD3KWz7FRmeIyEgRqSUiyV5iLZFXud5HRCJwiTGHMnyWpnKwBGEq0odyfHv9if4eqKpHcAnhXGAH8AxwraquPMkY3sBV/r7jfREWugbIFJF9wC248nlEpJkXazM/z/8mULTM/fe4X+O7cRW7bxQ9qJwuL/K5HhCR+j77PwAW4BLVx7gKeYCXgNeAGcB63Jf5rwC8+orzgF/jmvFmAN38iKUO8DzuWjfgKqgfK8/FmdARmzDImOpLXIe9Nqq6JtSxmKrH7iCMMcYUyxKEMcaYYlkRkzHGmGLZHYQxxphi/WxQtKqqXr16mpqaGuowjDGmSlmwYMEOVU0pbl+1SRCpqanMnz+/9BcaY4w5SkQ2nGifFTEZY4wpliUIY4wxxbIEYYwxpljVpg7CGBNcubm5ZGVlkZOTU/qLTaUTHR1NkyZNiIiI8PsYSxDGGL9kZWURHx9PamoqIhLqcMxJUFV27txJVlYWLVq08Ps4K2IyxvglJyeH5ORkSw5VkIiQnJx80nd/liCMMX6z5FB1leXfrsYniCN5BfzlkxVk7c4u/cXGGFOD1PgEsXVfDm98t5Hbxi8kJzc/1OEYY05g586dpKWlkZaWRsOGDWncuPHR9SNHjpR47Pz587nzzjtLfY9TTz01ILFOnz6dhISEo/GlpaXx5ZdfBuTcFanGV1I3TYrlHyO7cdNrC/jDh8v5yyVdQh2SMaYYycnJZGRkAPDII48QFxfHvffee3R/Xl4etWoV/5WWnp5Oenp6qe8xa9asgMQKMHDgQD766KMT7ldVVJWwsLBi10+kpOsMtBp/BwFwVqeG3DakFW/O3cjb80qb994YU1mMHj2aW265hT59+nD//fczd+5c+vXrR/fu3Tn11FNZtWoV4H7RX3DBBYBLLtdffz1DhgyhZcuWPPnkk0fPFxcXd/T1Q4YMYcSIEbRv356rrrqKwpGvp0yZQvv27enZsyd33nnn0fP6IzMzk3bt2nHttdfSuXNnZs6cedz6pk2buO++++jcuTNdunThrbfeOhrPwIEDGTZsGB07dgzIZ+ePGn8HUejXZ7VjSdZeHvzgezo2qkPnxgmhDsmYSusPHy5j+ZZ9AT1nx0Z1ePjCTid9XFZWFrNmzSI8PJx9+/Yxc+ZMatWqxZdffsn//d//8d577/3smJUrVzJt2jT2799Pu3btuPXWW3/WP2DRokUsW7aMRo0a0b9/f7799lvS09O5+eabmTFjBi1atGDUqFEnjGvmzJmkpaUdXX/vvfcIDw9n9erVvPrqq/Tt25fMzMzj1t977z0yMjJYvHgxO3bsoFevXgwaNAiAhQsX8v33359UM9XysjsIT3iY8O8r0qhXO5JbXl/AnuySyzSNMZXDZZddRnh4OAB79+7lsssuo3PnzowZM4Zly5YVe8z5559PVFQU9erVo379+mzduvVnr+nduzdNmjQhLCyMtLQ0MjMzWblyJS1btjz6JV1Sghg4cCAZGRlHl1atWgHQvHlz+vbte/R1vuvffPMNo0aNIjw8nAYNGjB48GDmzZt3NJ6KTA5gdxDHSY6L4pmrezLyv7O5a0IGL4/uRViYNeszpqiy/NIPltq1ax99/vvf/56hQ4cyceJEMjMzGTJkSLHHREVFHX0eHh5OXl5emV5T3niLW/f3uIpgdxBFpDVN5OFhHfn6h+08MXV1qMMxxpyEvXv30rhxYwBeeeWVgJ+/Xbt2rFu3jszMTICjdQSBMnDgQN566y3y8/PZvn07M2bMoHfv3gF9j5NhCaIYV/ZuxoieTXhy6mo+/f6nUIdjjPHT/fffz29/+1u6d+8esF/8vmJiYnjmmWc455xz6NmzJ/Hx8SQkFF9fWVgHUbi8++67pZ7/4osvpmvXrnTr1o3TTjuNv//97zRs2DDQl+G3ajMndXp6ugZywqCc3HyuGPcdq37az7u39qNTI6u0NjXbihUr6NChQ6jDCLkDBw4QFxeHqnL77bfTpk0bxowZE+qw/FLcv6GILFDVYtsA2x3ECURHhDPu2p4kxkZw46vz2bbfRrA0xsDzzz9PWloanTp1Yu/evdx8882hDiloLEGUoH58NM9fm87u7Fxu+t8C62ltjGHMmDFkZGSwfPlyxo8fT2xsbKhDChpLEKXo3DiBf13ejYxNe3jgvSVUlyI5Y4wpjSUIP5zT+RTuPastkzK28Mz0taEOxxhjKoT1g/DT7UNbs3rbAR77bBWtUuI4p3PoWhYYY0xFsDsIP4kIf7u0K2lNE7lzwiJmr90Z6pCMMSaoLEGchOiIcF4e3YvmSbH88tV5LN60J9QhGVNjDB06lM8+++y4bU888QS33nrrCY8ZMmQIhc3fzzvvPPbs2fOz1zzyyCM8/vjjJb73pEmTWL58+dH1hx56KCDDd1f2YcEtQZykurUjee2GPiTFRXLdy3P5Yev+UIdkTI0watQoJkyYcNy2CRMmlDgekq8pU6aQmJhYpvcumiAeffRRzjjjjDKdq6iiYzYVPa+qUlBQcML1EwlER0FLEGXQMCGa12/oQ0R4GNe8OIdNu2w2OmOCbcSIEXz88cdHJwfKzMxky5YtDBw4kFtvvZX09HQ6derEww8/XOzxqamp7NixA4CxY8fStm1bBgwYcHRIcHB9HHr16kW3bt249NJLyc7OZtasWUyePJn77ruPtLQ01q5dy+jRo4/2jJ46dSrdu3enS5cuXH/99Rw+fPjo+z388MP06NGDLl26sHLlSr+vtbIMC26V1GXUPLk2r9/Qh5HPzeaqF+bw7i39qF8nOtRhGVMx7r4bvMl7AiYtDZ544oS7k5KS6N27N5988gnDhw9nwoQJjBw5EhFh7NixJCUlkZ+fz+mnn86SJUvo2rVrsedZsGABEyZMICMjg7y8PHr06EHPnj0BuOSSS7jxxhsBePDBB3nxxRf51a9+xbBhw7jgggsYMWLEcefKyclh9OjRTJ06lbZt23Lttdfy7LPPcvfddwNQr149Fi5cyDPPPMPjjz/OCy+88LN4KvOw4EG9gxCRc0RklYisEZEHitl/j4gsF5ElIjJVRJp729NEZLaILPP2XR7MOMuqXcN4XvlFL3YcOMzVL85h90EbItyYYPItZvItXnr77bfp0aMH3bt3Z9myZccVBxU1c+ZMLr74YmJjY6lTpw7Dhg07uu/7779n4MCBdOnShfHjx59wuPBCq1atokWLFrRt2xaA6667jhkzZhzdf8kllwDQs2fPowP8FVWZhwUP2h2EiIQDTwNnAlnAPBGZrKq+/3KLgHRVzRaRW4G/A5cD2cC1qrpaRBoBC0TkM1XdE6x4y6p7s7q8cG06o1+Zx+iX5/L6L/sQHx1R+oHGVGUl/NIPpuHDhzNmzBgWLlxIdnY2PXv2ZP369Tz++OPMmzePunXrMnr0aHJyyjY0zujRo5k0aRLdunXjlVdeYfr06eWKt3DI8LIMF14ZhgUP5h1Eb2CNqq5T1SPABGC47wtUdZqqFhbgfwc08bb/oKqrvedbgG1AShBjLZdTW9fjmSt7sGzLPm54ZT6HjtiQHMYEQ1xcHEOHDuX6668/evewb98+ateuTUJCAlu3buWTTz4p8RyDBg1i0qRJHDp0iP379/Phhx8e3bd//35OOeUUcnNzGT9+/NHt8fHx7N//8wYp7dq1IzMzkzVr1gDw2muvMXjw4EBcaokqaljwYCaIxoDvBM9Z3rYTuQH42b+siPQGIoGfdWEWkZtEZL6IzN++fXs5wy2fMzo24F+XpzFvwy5ufn0Bh/MsSRgTDKNGjWLx4sVHE0S3bt3o3r077du358orr6R///4lHt+jRw8uv/xyunXrxrnnnkuvXr2O7vvjH/9Inz596N+/P+3btz+6/YorruCxxx6je/furF177KsoOjqal19+mcsuu4wuXboQFhbGLbfcclLXU5mHBQ/acN8iMgI4R1V/6a1fA/RR1TuKee3VwB3AYFU97LP9FGA6cJ2qflfS+wV6uO+yenveJu5/bwlnd2rA01f2oFa4NRQz1YMN9131VabhvjcDTX3Wm3jbjiMiZwC/A4YVSQ51gI+B35WWHCqTkb2a8vCFHfls2VbufWcxBQU2uJ8xpmoKZjPXeUAbEWmBSwxXAFf6vkBEugPP4e40tvlsjwQmAv9T1dLvtyqZX/RvQfaRfB77bBUxkeGMvaiLzW1tjKlygpYgVDVPRO4APgPCgZdUdZmIPArMV9XJwGNAHPCOiABsVNVhwEhgEJAsIqO9U45W1YxgxRtotw9tTfaRPJ6etpac3AL+PqIrEVbcZKo4VcX7v2qqmLJUJwS1o5yqTgGmFNn2kM/zYvuqq+rrwOvBjK0i3HtWO2Ija/HYZ6vYdyiXp6/qQXREeKjDMqZMoqOj2blzJ8nJyZYkqhhVZefOnURHn1xnXutJHUQiwu1DW5MQE8HvP/iea1+aywvXpVPH+kmYKqhJkyZkZWUR6haDpmyio6Np0qTJSR1jCaICXN23OQkxEYx5K4NR477j1et7Uy8uKtRhGXNSIiIiAtZD11QNViheQS7s1ojnr0tn7fYDjPzvbLJ22wB/xpjKzRJEBRrarj6v39CH7QcOc+mzs1i2ZW+oQzLGmBMqMUGISLiIjKmoYGqC9NQk3r3lVMJEuPy575i52spzjTGVU4kJQlXzAf9m4zB+a9cwnom39adJ3Rh+8fI83pm/qfSDjDGmgvlTxPStiDwlIgNFpEfhEvTIqrmGCdG8c0s/+rZM5r53l/Dk1NVlaqdsjDHB4k8rpjTv8VGfbQqcFvBoapj46AheGt2LB95fwj+/+IEtew7xx4s6W4c6Y0ylUGqCUNWhFRFITRVZK4x/XNaNxokx/OerNWTtPsTTV/UgIcb6ShhjQqvUn6oikiAi/ywcVltE/iEiCRURXE0hIvz6rHY8NqIrc9bv5JJnvmXDzoOhDssYU8P5U5bxErAfNz7SSGAf8HIwg6qpLktvyus39GHnwSNc9PS3zFm3M9QhGWNqMH8SRCtVfdibGW6dqv4BaBnswGqqPi2TmXRbf+rWjuTqF+dYCydjTMj4kyAOiciAwhUR6Q8cCl5IJrVebSbe2p/eLZK4790l/OWTFeTbvBLGmArmTyumW4D/+dQ77AauC15IBiAhNoJXftGbRyYv47mv17Hix/08eUUaibGRoQ7NGFNDlNqTGrhGVbsBXYGuqtpdVZdUSHQ1XER4GGMv7sJfLunCd2t3cuFT37B8y75Qh2WMqSH86Uk9wHu+T1Xt2ykERvVuxls39yU3T7nk2W/5IONnM7caY0zA+VMHsUhEJovINSJySeES9MjMcbo3q8uHvxpA18aJ3DUhgz9+tJy8/IJQh2WMqcb8SRDRwE5cz+kLveWCYAZlipcSH8X4G/sw+tRUXvxmPZeP+45Nu2zYcGNMcJRYSe3VQexU1XsrKB5TiojwMB4Z1onuzRJ5cOL3nPfvmYy9pAvDujUKdWjGmGrGnzqI/hUUizkJw9MaM+WugbRpEMedby7i3ncWc+BwXqjDMsZUI/40c80QkcnAO8DR8R9U9f2gRWX80jQplrdv7seTX63hqa9WMz9zF/++ojvdmiaGOjRjTDVgdRBVXK3wMO45sy0TburHkbwCLn12Fk9OXW0V2MaYcpPqMgdBenq6zp8/P9RhhNTe7Fwemvw9H2RsoVvTRP45shutUuJCHZYxphITkQWqml7cvhPeQYjI2z7P/1Zk3+eBC88ESkJsBP++ojtPXdmdzB0HOf/Jmbw6K5MCG6bDGFMGJRUxtfF5fmaRfSlBiMUEyAVdG/H5mEH0aZHMw5OXcd3Lc/lxrw2fZYw5OSUliJJ+dtpP0kquQZ1oXvlFL/50UWfmZ+7mrH/O4LXvNtjdhDHGbyUliFgR6S4iPYEY73mPwvUKis+Ug4hwdd/mfHr3QLo2TeD3k77nsudm88PW/aEOzRhTBZywklpEppV0YGWbitQqqUumqry/cDN/+ng5Bw7ncevgVtw2tDXREeGhDs0YE0IlVVJbK6YaZueBw4z9eAXvL9pMy3q1+dPFnTm1Vb1Qh2WMCZEytWIy1VNyXBT/vDyN127oTV6BcuXzc7h7wiK27c8JdWjGmErGEkQNNbBNCp+PGcSdp7VmytKfOP3xr3l1VqbNXGeMOcoSRA0WHRHOPWe149O7B9KtaSIPT17G8Ke/IWPTnlCHZoypBEpNEOJcLSIPeevNRKS3PycXkXNEZJWIrBGRB4rZf4+ILBeRJSIyVUSa++y7TkRWe4tNcRpELVPieO2G3jx1ZXe27TvMxc98y2/eXcLOA4dDHZoxJoRKraQWkWeBAuA0Ve0gInWBz1W1VynHhQM/4DrZZQHzgFGqutznNUOBOaqaLSK3AkNU9XIRSQLmA+m4PhcLgJ6quvtE72eV1IGxPyeXJ6eu5uVvM4mJDOeeM9tyTd/m1Aq3m01jqqPyVlL3UdXbgRwA70s60o/jegNrVHWdqh4BJgDDfV+gqtNUtXDGm++AJt7zs4EvVHWX935fAOf48Z6mnOKjI/jd+R359O6BpDVN5A8fLuf8J79h9tqdoQ7NGFPB/EkQud7dgAKISArujqI0jYFNPutZ3rYTuQH45GSOFZGbRGS+iMzfvn27HyEZf7WuH8//ru/Nf6/uycEjeYx6/jtufX0Ba7ZZJztjagp/EsSTwESgvoiMBb4B/hzIIETkalxx0mMnc5yqjlPVdFVNT0mx4aECTUQ4p3NDvrxnMGPOaMuMH7Zz1r9mcO87i8nabVOdGlPdlTblaBiwHrgfOB0Q4CJVXeHHuTcDTX3Wm3jbir7HGcDvgMGqetjn2CFFjp3ux3uaIIiOCOeuM9pwdd9mPDt9Lf/7bgMfZGzmqj7NuX1oa1Lio0IdojEmCPyppF6kqt1P+sQitXCV1KfjvvDnAVeq6jKf13QH3gXOUdXVPtuTcBXTPbxNC3GV1LtO9H5WSV1xftx7iCenrubt+VlEhodxbb/m3DioJfXiLFEYU9WUt5J6qohcKiJyMm+qqnnAHcBnwArgbVVdJiKPisgw72WPAXHAOyJSOLUpXiL4Iy6pzAMeLSk5mIp1SkIMf7mkK1/eM5izOjXg+ZnrGPi3afx5ygp2WNNYY6oNf+4g9gO1gTxcSyYBVFXrBD88/9kdROis3X6Ap75awwcZm4msFcbVfZpz8+BWVvRkTBVgg/WZCrHOSxSTMjYTER7GqN7NuHlwS05JsNHhjamsyp0gvM5xbYDowm2qOiNgEQaAJYjKY/2OgzwzbQ0TF20mTIQR6U24dXArmibFhjo0Y0wR5UoQIvJL4C5cS6IMoC8wW1VPC3Cc5WIJovLZtCub52as5e15WeSrclFaY24d0orW9eNCHZoxxlPeBLEU6AV8p6ppItIe+LOqXhL4UMvOEkTltXVfDuNmrGP8nA3k5BZwRof63DSoFb1S63KSbR+MMQFW3lZMOaqa450oSlVXAu0CGaCp3hrUieb3F3Tk29+cxl2nt2HBht2MfG42Fz8zi0+W/mhDjBtTSZXYUc6TJSKJwCTgCxHZDWwIZlCmekqOi2LMmW25ZXAr3l2wiednrufW8QtJTY7lulNTGdGzCfHREaEO0xjjOalWTCIyGEgAPvUG4Ks0rIip6skvUD5b9hPPz1zHoo17qB0ZzoieTbj21FRapVg9hTEVobx1EM2K266qGwMQW8BYgqjalmTt4ZVZmXy0+EeO5BcwsE09Rp+aypB29QkPs3oKY4IlEJXUiusgFw20AFapaqdAB1oeliCqhx0HDvPmnI28PmcDW/cdpmlSDFf1ac7I9KYk1fZnlHljzMkIaEc5EekB3KaqvwxEcIFiCaJ6yc0v4PNlW/nf7EzmrN9FZK0wLuh6Ctf2S6VbkwRr/WRMgAS8J7WILFXVLuWOLIDKnCB27IDRo+GRRyC92M/IhNgPW/fz2uwNvL8wi4NH8mnfMJ6R6U25qHtju6swppzKW8R0j89qGG6E1WRVPTtwIZZfuRJEz57u+YIFUK9eYAMzAXPgcB4TF23mnfmbWJK1l4hw4fT2DRjZqwmD2qTYtKjGlEF5E8TDPqt5QCbwXmHfiMqiXEVM8+bBgAEwaBB8+imEhwc2OBNwK3/axzvzs5i0aDM7Dx4hJT6K4d0acXGPxnQ8pY4VQRnjJxuszx8vvAA33gi//S38OaAT5pkgOpJXwLRV23h3QRbTV20jN19p1yCei7o35qLujWygQGNKUd47iA/x5qMujqoOO9G+ihSQSuobb3SJYuJEuOiigMRlKs7ug0f4aOmPTFyYxcKNexCBPi2SuLBbI87tfIrVVxhTjPImiH8DDYHXvU2jgK24ntWo6tcBi7QcApIgcnJg4ED44QdX7NS2bWCCMxVuw86DTFy0mQ8Xb2Ht9oOEhwmntkrmwq6NOLtTQxJirce2MVD+BDG/6MHFbQu1gDVz3bDBVVo3aABz5kCc9eitylSVFT/u56MlW/hoyY9s3JVNRLgwuG0Kw9Mac0aHBsREWp2TqbnKmyBWAOer6jpvvQUwRVU7BDzScghoP4gvv4Szz4YRI2DCBLAKz2pBVVmStZcPF2/hwyVb2LrvMLUjwzm7c0OGpzWmf6tkawllapzyJohzgHHAOlxv6ubATar6eaADLY+Ad5T7619dhfVf/gIPPBC485pKIb9AmbNuJx9kbGHK9z+yPyeP5NqRDG6XwtB29RnUJsWKoUyNEIgZ5aKA9t7qSlWtdDPTBzxBqMKVV8Jbb7lK6+HDA3duU6nk5OYzfdU2piz9iRmrt7MnO5fwMKFns7oMaZ/C4LYp1nTWVFtlShAi0gvYpKo/eevXApfihvp+RFV3BSneMgnKUBuHDrm+EStWwKxZ0LVrYM9vKp38AiVj0x6mrdzGtFXbWLZlHwD14qIY1KYeg9qmMKBNPerFRYU4UmMCo6wJYiFwhqruEpFBwATgV0Aa0EFVRwQp3jIJ2lhMmzdDr14QGelaNqWkBP49TKW1dV8OM37YzozVO/hm9XZ2Z+cC0LlxHc7q2JBzOzekTYP4EEdpTNmVNUEsVtVu3vOnge2q+oi3nqGqacEJt2yCOljf3LkweLBLFF9+6ZKFqXHyC5RlW/Yy44ftTFu1nYUbd6MKrevHcV7nhpzb5RTaN4y3oihTpZQ1QXwPpKlqnoisxFVMzyjcp6qdgxZxGQR9NNc333R1EjfcAM8/by2bDFv35fDZsp+YsvRH5q7fRYFCk7ox9G2ZTO/UJHq3SKJ5cqwlDFOplZQgSppy9E3gaxHZARwCZnonaw3sDXiUld2oUbBsGYwdC506wZgxoY7IhFiDOtFc2y+Va/ulsuPAYT5ftpVpq7YxdcVW3l2QBUBKfBS9U5Po2yqZAa3rkWoJw1QhJbZiEpG+wCnA56p60NvWFohT1YUVE6J/KmQ+iIICuOwyeP99GD/e3VEYU0RBgbJ2+wHmZu5i7nq3/LjXjW3ZODGG/q2T6d+6Hv1bW2W3Cb2ADdYnIjep6riARRZAFTZhUE4OnHMOfPstTJ4M554b/Pc0VZqqkrkzm2/W7ODb1TuYtXYH+3LyAGjfMJ5+rZI5tVU9erdIIiHG+l6YihXIBLFQVXsELLIAqtAZ5fbuhSFDYNUqmDoV+vWrmPc11UJ+gfL95r18s2YH363bybzMXeTkFhAm0KVxAn1buTqM9OZJ1lnPBF0gE8QiVe0esMgCqMKnHN261c0hsXMnzJzp6iWMKYPDeflkbNzDrLU7mbV2Bxmb9pCbr4hAuwbx9G6RRC+v0rtBnehQh2uqmUAmiCaqmhWwyAIoJHNSr18P/ftDWJgrcmrevGLf31RLObn5ZGzaw7z1u5ibuYuFG3Zz8Eg+AM2SYundIoneqUn0apFkld6m3Mo7FlMUrgd1Kj6tnlT10QDGWG4hSRAAS5e63tb167s7ifr1Kz4GU63l5Rew/Md9zF2/i3mZu5iXuZtdB48Arod3z+aJ9Gxelx7N6tK5cQLRETY6rfFfeRPEp7hmrQuA/MLtqvqPQAZZXiFLEODuHs48E1q2hGnTrLe1CSpVr5XU+t3My9zFwo272bAzG4CIcKFTowR6Nq9LevO69EytS/14K5YyJ1beBFHmTnHeSLD/BsKBF1T1r0X2DwKeALoCV6jquz77/g6cD4QBXwB3aQnBhjRBAHz1FVxwAbRu7Z7Xqxe6WEyNs33/YRZu3O2WDbtZnLWXI3kFADRPjvUSRhK9UuvSKiWOsDArljJOWTvKFZolIl1UdelJvmk48DRwJpAFzBORyaq63OdlG4HRwL1Fjj0V6I9LHADfAIOB6ScTQ4U67TT48EOXJE4/3SWJ5ORQR2VqiJT4KM7u1JCzOzUEXMX395v3sWDDLuZn7ubrVdt5f+FmABJiIo7eXfRKTaKLFUuZE/AnQQwARovIeuAwbk4IVdXShjbtDazxmWhoAjAcOJogVDXT21dQ5FgFooFI7/0icNOcVm6nn+76Rlx4IZxxhhu3yZKECYGoWuH0bF6Xns3rctOgY30x5me6hDF/wy6mrtwGuGKpDqfUoWuTBLo1SaRb00RapcQRbncZNZ4/CaKsPcEaA5t81rOAPv4cqKqzRWQa8CMuQTylqiuKvk5EbgJuAmjWrFkZwwywM8+EDz5w80eceaZLEklJoY7K1HAiQot6tWlRrzaXpTcFYNfBIyzYsJsFG3azJGsPkxZt4fXvNgJQOzKcNg3iaZ4cS/OkWJomxdI8uTapybGkxEdZy6kaotQEoaobAESkPu5XfdB54z11AJp4m74QkYGqOrNIbONws92Rnp7uf3vdYDv7bJg0ySWJs85ySSIxMdRRGXOcpNqRnNmxAWd2bAC4IULW7TjI4k17WJK1hzXbD7Bgw24+XLyFAp//XSnxUXRvmkhas0TSmibStUkicVH+/NY0VU2p/6oiMgz4B9AI2IabcnQFUFrPsM1AU5/1Jt42f1wMfKeqB7wYPgH64Q0YWCWcc46bie6ii9zzL76AeJs3wFReYWFC6/pxtK4fx6U9mxzdnptfwObdh9iwK5v12w+wJGsvizbt4fPlrtQ3TNyQ5+0b1qH9KfG0bxhPu4Z1aJQQbXcaVZw/af+PQF/gS1XtLiJDgav9OG4e0EZEWuASwxWAv6PbbQRuFJG/4IqYBuNaO1Ut550Hb78NI0bA+efDJ59A7dqhjsqYkxIRHkZqvdqk1qvN4LbHmnDvPniEjKw9ZGzcw9LNe1mwYTeTF285ur9OdC1a1Y+jRbI7tnlyLC2889SJtiFEqgJ/mrnOV9V0EVkMdFfVAt/JhEo59jzcF3s48JKqjhWRR4H5qjrZm9Z0IlAXyAF+UtVOXguoZ4BBuArrT1X1npLeK+TNXEvy1ltu5NchQ+CjjyAmJtQRGRMUew/l8sPW/az8aT+rftrHuu0HydxxkC3eaLaFGtSJonX9ONrUj6dV/Thap8TRol5t6sVFUis8LETR10zl7QfxJXAR8FcgGVfM1EtVTw1wnOVSqRMEwGuvwXXXHSt6irJhnk3NkZObz4ad2azfcZD1Ow6yZtsB1mw/wJqt+48OIwJuHq7k2pGkxEdTPz6K+vFRNK4bQ7Ok2KOLVZIHVnn7QQzHTRh0N3AVkABUqmE2qoRrrnFDhd90E1x+ObzzDkTYbbapGaIjwmnXMJ52DY+vh1NVftqXw+qtB9i4K5vt+w+zbf9htu/PYdv+w6z8aR/b9h/G93dsVK0wGteNISk2koSYCBJiI0iMiSQxNoJ6cVE0TYqhSd1YGifGEFnL7kbKw59WTAdFpDnQRlVfFZFYXJGROVk33ghHjsAdd7gk8eabdidhajQR4ZSEGE5JOHGx6+G8fDbvPsTGXdls2pXNhp3ZbNl7iD3Zufy4N4eVP+1n76FcDhzOK3JuaFgnmiZ1Y0iuHUXd2pHUjY2gbqxLJomxkcRH1yIuqhZ1oiOI855bUjnGn1ZMN+L6GiQBrXD9G/4LnB7c0Kqp22+H/Hy46y7XDPb99yE2NtRRGVNpRdUKp2VKHC1T4kp8XW5+Adv2HyZrVzabdh8ia3c2m3a5x7XbD7B7Qy57so+QV1BysXqtMCEiPIyIcCGyVhgR4WFE1gqjbmwkKfFRbomLol58FIkxEdQKE0SE8DAhPAzCRAgTOTptveCei0CtsDDCw4SIcPf6WmHufaIiwokMDyMqIoyoWmFEevUwBermDylQPfro3kcI994zmMVt/hQx3Y7rFT0HQFVXe30iTFndeSfExcEvf+nqJD76COrUCXVUxlRpEeFhNE6MoXFizAl75KoqBw7nsftgLnsOHeFATh77D+exPyePAzm57M/JIycvn9x85UheAbn5bsnJLWDXwSNs2pXNwg272ZV9hJOYKSGowsOE7k0TeffWwFcL+5MgDqvqkcIsJSK1cC2LTHlcf71r8nr11W6Ijk8/tWE5jAkyESE+OoL46AiaUfY799x8lzD2Hso9+su+oADyvV/6hY1/FFB1ialAoUCV3PwC8guU3Hz1Hgs4klfA4bx8DucVHF0A7y7B9VEJ9+5MCt8jv0DJK1AKCpQGCcHpw+xPgvhaRP4PiBGRM4HbgA+DEk1Nc/nlLkmMGOGawH7xBTRsGOqojDGliAgPo0Gd6Go/w58/tTEPANuBpcDNwBTgwWAGVaNccAFMmeJmpxs4ENatC3VExhgD+JEgVLVAVZ9X1ctUdYT33IqYAum009x4Tbt2Qd++MHt2qCMyxpgTJwgRWVLSUpFB1giFiSEhAYYOdUN0GGNMCJVUB1GAq2N5A1fncKhCIqrJ2rZ1SeKii1z9xLp18JvfgPUaNcaEwAnvIFQ1DRgFxOGSxFjcCK6bC4cAN0FQr54rbho1Cn77W9e5Ljc31FEZY2qgEusgVHWlqj6sqj1wdxH/A8ZUSGQ1WXQ0jB8Pv/89vPginHsu7N4d6qiMMTVMic1cRaQxbpjui4HduOQwsQLiMiLw6KPQqpW7i+jXz3Woa9061JEZY2qIkiqpv8bdNUQAvwCuAz4GIkXE5tCsKNddB1Onwo4d0KcPTJ8e6oiMMTVESUVMzXHzNNwMfAbM95YF3qOpKAMHwpw50KCBm+f6xRdDHZExpgY4YRGTqqZWYBymNK1auRZOI0e6MZxWroS//hXCbWBdY0xw2Li2VUlCAnz8sRsR9vHH3ZSm27aFOipjTDVlCaKqqVULnnoKxo2Dr7+GtDT3aIwxAWYJoqq68UZXLxEf74bqGDsWCgpCHZUxphrxK0GIyAAR+YX3PEVEWgQ3LOOXbt1g/nzX6/rBB93cElu3hjoqY0w1UWqCEJGHgd8Av/U2RQCvBzMocxLi412nunHjYOZMV+T0zTehjsoYUw34cwdxMTAMOAigqluA+BKPMBVL5FiRU1ycG+zvmWeoNFNeGWOqJH8SxBFveG8FEJHawQ3JlFnXrjBvHpx9tmvpdP31kJMT6qiMMVWUPwnibRF5DkgUkRuBL4HngxuWKbPERJg8GR56CF55BQYMgI0bQx2VMaYK8mfCoMeBd4H3gHbAQ6r6n2AHZsohLAz+8Af44ANYvRp69nTTmRpjzEnwqxWTqn6hqvep6r2qat80VcWwYTB3LqSkwFlnwZVXwpYtoY7KGFNF+NOKab+I7CuybBKRiSLSsiKCNOXQrh0sWOCKnN5/363/8582x4QxplT+3EE8AdwHNAaaAPfiJhCaALwUtMhM4MTEuCKnZctg8GD49a+he3cbGdYYUyJ/EsQwVX1OVfer6j5VHQecrapv4UZ7NVVFq1ZuTonJk+HgQdcc9oYbYN++UEdmjKmE/EkQ2SIyUkTCvGUkUNh20hraV0UXXgjLl7spTV95Bbp0cXNOGGOMD38SxFXANcA2YKv3/GoRiQHuCGJsJphiYuDPf4ZZs9zzM86AO+5wdxbGGIN/zVzXqeqFqlpPVVO852tU9ZCq2pgOVV2fPrBoEYwZ43pfd+sG334b6qiMMZWAP62YokXkdhF5RkReKlz8ObmInCMiq0RkjYg8UMz+QSKyUETyRGREkX3NRORzEVkhIstFJNXvqzInJybGtWyaNs2NCDtwINx5J+zfH+rIjDEh5E8R02tAQ+Bs4GtcS6ZSvzlEJBx4GjgX6AiMEpGORV62ERiNaxVV1P+Ax1S1A9AbV8RlgmnwYFi82BU1PfUUdOoEU6aEOipjTIj4kyBaq+rvgYOq+ipwPtDHj+N6A2u8IqojuGaxw31foKqZqroEOG4iAy+R1CrslKeqB1Q124/3NOUVHw9PPumKmeLj4fzzXQe77dtDHZkxpoL5kyAKe1TtEZHOQAJQ34/jGgObfNazvG3+aOu93/siskhEHvPuSI4jIjeJyHwRmb/dvsACq18/WLgQHnkE3n0XOnSAF1+E/PxQR2aMqSD+JIhxIlIXeBCYDCwH/hbUqKAWMBDXKa8X0BJXFHUcVR2nqumqmp6SkhLkkGqgqCh4+GHIyHA9sH/5SzffxJQpNpS4MTVAiQlCRMKAfaq6W1VnqGpLVa2vqs/5ce7NQFOf9SbeNn9kARle8VQeMAno4eexJtA6dnSTEL3zDhw65IqdTj/dzWZnjKm2SkwQqloA3F/Gc88D2ohICxGJBK7A3YH4e2yiiBTeFpyGu3MxoSICI0a4Dnb/+Q8sXQq9ern6CRtO3JhqyZ8ipi9F5F4RaSoiSYVLaQd5v/zvAD4DVgBvq+oyEXlURIYBiEgvEckCLgOeE5Fl3rH5uOKlqSKyFBBsDorKITLStXJauxZ+9zuYNAnat4exY21yImOqGdFSypJFZH0xm1VVK9VIrunp6Trfijwq3saNcM898N57bqynf//bFUEZY6oEEVmgqunF7fOnJ3WLYpZKlRxMCDVr5lo5ff45RETABRe4sZ7Wrg11ZMaYcvKnJ3WsiDwoIuO89TYickHwQzNVyplnuk52jz3mhhHv0AFuuw02+9suwRhT2fhTB/EycAQ41VvfDPwpaBGZqisyEu69F1atcsOIP/+8K3YaMwa2bg11dMaYk+RPgmilqn/H6zDn9WiWoEZlqrZGjeDZZ+GHH1wrp//8B1q2hN/8BnbuDHV0xhg/+ZMgjnhDeyuAiLQCDgc1KlM9tGgBL73kmsZedJErfmrVCh5/3Fo8GVMF+JMgHgE+BZqKyHhgKmXvG2FqorZtYfx4WLIE+veH++5zTWPffNONHmuMqZT8acX0OXAJbqiLN4F0VZ0e3LBMtdS5M3z8MXzxBSQmuuKnvn1hxoxQR2aMKYY/rZg+BM4CpqvqR6q6I/hhmWrtjDNgwQI33emWLW6Y8VNPhTfegCNHQh2dMcbjTxHT47iB85aLyLsiMkJEooMcl6nuwsPhuutcRfYTT8COHXDVVa5fxUMPWfNYYyoBf4qYvlbV23Ajqj4HjMQm7zGBEhsLd90FK1fCZ59B797wpz9B8+Zw+eXuTsMYExL+3EHgtWK6FLgFN/z2q8EMytRAYWFw1lkwebLrhT1mjEsY6elu5NjPP7chxo2pYP7UQbyNG2zvNOApXL+IXwU7MFODtWjhmsRu3OgeV66Es8+GHj1cy6e8vFBHaEyN4M8dxIu4pHCLqk4DThWRp4MclzFQp47rmb1unetPkZPjWj61aAF//rNNg2pMkPlTB/EZ0FVE/i4imcAfgZXBDsyYo6Ki4Be/gGXLXBFU+/ZuqPGmTWH0aKunMCZITpggRKStiDwsIiuB/+DmlxZVHaqq/6mwCI0pFBbmRor94gvXO/uGG9xIsunpbg7t//7XhvIwJoBKuoNYiat3uEBVB3hJwWasN5VDhw7w9NOuOewTT8DevXDrrdCwoUsiEyZAdnaoozSmSispQVwC/AhME5HnReR0bJA+U9kkJLhmssuWwaJFcPfd7nHUKGjQwBVBffuttYAypgxOmCBUdZKqXgG0B6YBdwP1ReRZETmrguIzxj8ikJbmWj1t2ABffeX6Ubz3HgwY4Ib5eOIJK4Iy5iT4U0l9UFXfUNULgSbAIuA3QY/MmLIKD4ehQ+GFF+DHH91jfLzrW9G4seux/emnkJsb6kiNqdRKnZO6qrA5qU2plixxkxi9/jrs2QPJyXDppe5OY/Bgl1iMqWHKNSe1MdVG165u8qKffoIPPnCd78aPdz21GzeGO+5w9RU2BLkxgCUIUxNFRcGwYS45bNsG77wDAwfCiy+6+ooWLeD++2HhQqvcNjWaJQhTs8XGwogRLkls2+aKn7p0gX/9C3r2dJ3yHnrIFU9ZsjA1jCUIYwrFx7sK7I8+csVQ48a5oqexY6FbNzcz3gMPwPz5lixMjWAJwpjiJCfDjTe65rI//gjPPQctW7r5tHv1OlYMlZFhycJUW9aKyZiTsWuXGw/qnXfcEOR5ea5X96hRbmndOtQRGnNSSmrFZAnCmLLascN1xHvjjWPzavfoAeee61pI9e0LERGhjdGYUliCMCbYsrLc+E8ffACzZ0N+vhuu/LTT4JxzXD+Ltm3dgIPGVCKWIIypSHv2uLqLTz91s+Jt3Oi216njRp7t1ctNrdq7NzRpEtJQjbEEYUyoqMKqVe6uYu5cmDfPNZktHOajQwc4/3w47zzo3x8iI0Mbr6lxLEEYU5nk5MDixa7X9iefwNdfu4QRHw9nnumGK7/4YjdSrTFBZgnCmMrswAGYOhWmTHFLVpa7kzj/fNcy6oILICYm1FGaaipkYzGJyDkiskpE1ojIA8XsHyQiC0UkT0RGFLO/johkichTwYzTmJCKi4Phw11fi40b4bvv3ORHs2fDyJFQvz5cc42bl3vuXDh4MNQRmxoiaHcQIhIO/ACcCWQB84BRqrrc5zWpQB3gXmCyqr5b5Bz/BlKAXap6R0nvZ3cQptrJz4fp0+HNN11z2j173HYR12mvc2fXw3vwYDflqt1lmDIo6Q6iVhDftzewRlXXeUFMAIYDRxOEqmZ6+342fKaI9AQaAJ8CxQZvTLUWHu5Gmj39dHd3sX49LF0K339/7PGjj+DRR90AhP36uWa1Q4e6FlJW4W3KKZgJojGwyWc9C+jjz4EiEgb8A7gaOKOE190E3ATQrFmzMgdqTKUXHu56abdu7SqwC+3fDzNnuma1X30FDz/sBheMinJ3F+npx5rWduhgc16YkxLMBFEetwFTVDVL5MTTYKvqOGAcuCKmCorNmMojPt41kT3vPLe+c6drFTV7thtU8LXX4Jln3L7YWJco+vVzS9++rn7DmBMIZoLYDDT1WW/ibfNHP2CgiNwGxAGRInJAVX9W0W2M8ZGcDJdc4hZwkx/98INLFvPmuQrwxx93Y0iBq8vo188Nbd6zJ3Tv7pKOMQQ3QcwD2ohIC1xiuAK40p8DVfWqwuciMhpIt+RgTBmEhbk5Ldq3h6uvdtsOHYIFC9xdxuzZMG2amzwJXAV427auWKpjR0hNPbY0bGhDhdQwQUsQqponIncAnwHhwEuqukxEHgXmq+pkEekFTATqAheKyB9UtVOwYjLG4Fo7DRjglkI//eSSRuEyffqxpFEoMhKaN3cDEvbt65bu3V19h6mWrKOcMaZ4Bw+6fhmZmW7ZsAHWrHF9MTZ57U8iI12S6NnTzZGRmuqSSGoq1Kvn7khMpRaqZq7GmKqsdm3X8qlDh5/v27wZ5sxxy3ffubuNvXuPf01sLJxyCqSkuKVePffYoAF06gRpae65qbTsDsIYExh79ri7jMK7jcxMV3S1ffvxS+FAheDqNbp1c8miUyd3F9KihUssVt9RIewOwhgTfImJbunW7cSvUXVNcZcudQMWZmS45Z//PD5xFNZ3pKa6u46EhOOXunVdEmncGBo1cncrJuAsQRhjKo6IK2oaOtQthY4cgXXr3F3H+vXHP65d64qv9u491jy3qMRElygaNHDJIynJPRYucXGuyCw29tgSFeWSUuFy5Ih7rFsXWrVy56jhLEEYY0IvMvJYc9wTUXVNdPfuhd27YcsWt2zefGzZvh1WrnRzh+/eDYcPlz2mwkTRurUr9oqNhVq13BIe7h4jItz2mJhjjzExbkj3XbuOLTt3ul7viYnH18cU3h1FRPx8iY527xFCliCMMVWDyLFf/6ec4vpplObQIfcFnZ19/HLwoEsehV/GkZHHnu/c6VprrVnj7l7mzoV33nGDJ5ZVVJTrgLh37/FFaaWpVev4BBQV5RKlqusEWVDgnnfrBhMnlj2+E719wM9ojDGVRUyMq6cor4IClyDy8o5/PHLEJaHs7GOP2dnufZOS3JKcfGykXVXYt8/d6ezY4R737Su+qCsn5+fnPnzYJcqwMLcUPm/duvzXWAxLEMYYU5rCL+SIiPKdR+RYRXuQvtQDydqRGWOMKZYlCGOMMcWyBGGMMaZYliCMMcYUyxKEMcaYYlmCMMYYUyxLEMYYY4plCcIYY0yxqs1w3yKyHdhQjlPUA3YEKJzKyK6v6qvu12jXFxrNVTWluB3VJkGUl4jMP9GY6NWBXV/VV92v0a6v8rEiJmOMMcWyBGGMMaZYliCOGRfqAILMrq/qq+7XaNdXyVgdhDHGmGLZHYQxxphiWYIwxhhTrBqfIETkHBFZJSJrROSBUMcTCCLykohsE5HvfbYlicgXIrLae6wbyhjLQ0Saisg0EVkuIstE5C5ve7W4RhGJFpG5IrLYu74/eNtbiMgc72/1LRGJDHWs5SEi4SKySEQ+8tar2/VlishSEckQkfnetir1N1qjE4SIhANPA+cCHYFRIuLHRLeV3ivAOUW2PQBMVdU2wFRvvarKA36tqh2BvsDt3r9bdbnGw8BpqtoNSAPOEZG+wN+Af6lqa2A3cEPoQgyIu4AVPuvV7foAhqpqmk//hyr1N1qjEwTQG1ijqutU9QgwARge4pjKTVVnALuKbB4OvOo9fxW4qCJjCiRV/VFVF3rP9+O+ZBpTTa5RnQPeaoS3KHAa8K63vcpeH4CINAHOB17w1oVqdH0lqFJ/ozU9QTQGNvmsZ3nbqqMGqvqj9/wnoEEogwkUEUkFugNzqEbX6BW/ZADbgC+AtcAeVc3zXlLV/1afAO4HCrz1ZKrX9YFL6p+LyAIRucnbVqX+RmuFOgBT8VRVRaTKt28WkTjgPeBuVd3nfoQ6Vf0aVTUfSBORRGAi0D60EQWOiFwAbFPVBSIyJMThBNMAVd0sIvWBL0Rkpe/OqvA3WtPvIDYDTX3Wm3jbqqOtInIKgPe4LcTxlIuIROCSw3hVfd/bXK2uEUBV9wDTgH5AoogU/qiryn+r/YFhIpKJK9Y9Dfg31ef6AFDVzd7jNlyS700V+xut6QliHtDGaz0RCVwBTA5xTMEyGbjOe34d8EEIYykXr7z6RWCFqv7TZ1e1uEYRSfHuHBCRGOBMXD3LNGCE97Iqe32q+ltVbaKqqbj/c1+p6lVUk+sDEJHaIhJf+Bw4C/ieKvY3WuN7UovIebjy0HDgJVUdG9qIyk9E3gSG4IYX3go8DEwC3gaa4YZFH6mqRSuyqwQRGQDMBJZyrAz7/3D1EFX+GkWkK64CMxz3I+5tVX1URFrifnEnAYuAq1X1cOgiLT+viOleVb2gOl2fdy0TvdVawBuqOlZEkqlCf6M1PkEYY4wpXk0vYjLGGHMCliCMMcYUyxKEMcaYYlmCMMYYUyxLEMYYY4plCcJUOSKiIvIPn/V7ReSRAJ37FREZUfory/0+l4nIChGZVmR7qogc8kYALVyuDeD7DikcPdWY0thQG6YqOgxcIiJ/UdUdoQ6mkIjU8hlLqDQ3ADeq6jfF7FurqmmBi8yYsrE7CFMV5eHm9x1TdEfROwAROeA9DhGRr0XkAxFZJyJ/FZGrvHkXlopIK5/TnCEi80XkB2/coMLB8x4TkXkiskREbvY570wRmQwsLyaeUd75vxeRv3nbHgIGAC+KyGP+XrSIHBCRf4mbI2KqiKR429NE5DsvromFcwyISGsR+VLcvBILfa4xTkTeFZGVIjLe65mO95ks987zuL9xmWpMVW2xpUotwAGgDpAJJAD3Ao94+14BRvi+1nscAuwBTgGicOP8/MHbdxfwhM/xn+J+PLXBjSoaDdwEPOi9JgqYD7TwznsQaFFMnI2AjUAK7m79K+Aib990IL2YY1KBQ0CGzzLQ26fAVd7zh4CnvOdLgMHe80d9rmUOcLH3PBqI9eLdixvrKAyYjUtWycAqjnWeTQz1v7MtoV/sDsJUSaq6D/gfcOdJHDZP3VwSh3HDZ3/ubV+K+2Iu9LaqFqjqamAdbiTVs4BrvSG45+C+UNt4r5+rquuLeb9ewHRV3a6u6Gk8MMiPONeqm2SmcJnpbS8A3vKevw4MEJEE3Jf51972V4FB3jhAjVV1IoCq5qhqtk+8WapagEtAqbikkYO7q7kEKHytqcEsQZiq7AlcWX5tn215eH/XIhIG+E5b6TuuT4HPegHH18cVHX9GAQF+5fOl3UJVCxPMwfJcRDmUdZwc388hHyisO+mNm7DnAtxdlKnhLEGYKkvdIGdvc/zUlJlAT+/5MNxsbCfrMhEJ88rsW+KKXj4DbvWGGUdE2nqjdJZkLjBYROqJm952FPB1KceUJIxjo51eCXyjqnuB3SIy0Nt+DfC1upn2skTkIi/eKBGJPdGJxc2tkaCqU3B1O93KEaepJqwVk6nq/gHc4bP+PPCBiCzG/Qouy6/7jbgv9zrALaqaIyIv4IpiFnqVutspZbpIVf1RRB7ADWMtwMeq6s/wzq28oqxCL6nqk7hr6S0iD+LmEbjc238d8F8vAawDfuFtvwZ4TkQeBXKBy0p4z3jc5xbtxXqPH3Gaas5GczWmihCRA6oaF+o4TM1hRUzGGGOKZXcQxhhjimV3EMYYY4plCcIYY0yxLEEYY4wpliUIY4wxxbIEYYwxplj/D0/Co/Ohi5ShAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "\n",
    "hist = pd.DataFrame(hist)#\n",
    "hist.rename(columns={0:\"Training Error\"},inplace=True)\n",
    "hist2 = pd.DataFrame(hist2)\n",
    "hist2.rename(columns={0:\"Validation Error\"},inplace = True)\n",
    "\n",
    "fig = sns.lineplot(data = hist)\n",
    "fig = sns.lineplot(data = hist2,palette=['r'])\n",
    "fig.set(xlabel=\"Number of Epochs\",ylabel=\"Average Mean-Square Error\",title = \"Error vs. Num Epochs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "574cf54f-5edc-421e-977f-e02528ac1dcf",
   "metadata": {},
   "source": [
    "After optimizing the parameters for this new neural network on daily data I am going to test this model on the test set now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1297,
   "id": "ae7267b9-edd5-4c67-b5b8-0ebd39b68118",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The model has loss of 0.3057524263858795 over the TEST SET\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(f\"\\nThe model has loss of {loss(model(x_test),y_test).item()} over the TEST SET\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1223,
   "id": "31a73964-1312-45b7-906b-e514ab996fdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The loss on the TEST SET from using the factors at time t as a proxy for the factors at time t+1 gives us an average MSE of 0.5822740197181702\n",
      "\n"
     ]
    }
   ],
   "source": [
    "lastVal = x_test[:,-1,:]\n",
    "lastSeqPredictionLoss = loss(lastVal,y_test).item()\n",
    "print(f'\\nThe loss on the TEST SET from using the factors at time t as a proxy for the factors at time t+1 gives us an average MSE of {lastSeqPredictionLoss}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1234,
   "id": "e9ddc133-5659-4c30-b033-9c62568953a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The loss on the TEST SET from using a moving average over the last 5 periods as a proxy for the factors at time t+1 gives us an average MSE of 0.31237014009948866\n",
      "\n"
     ]
    }
   ],
   "source": [
    "temp_X_Test = np.array(x_test[:,0,:])\n",
    "MALength = 30\n",
    "MA = np.zeros((temp_X_Test.shape[0]-MALength,5))\n",
    "\n",
    "for i in range(temp_X_Test.shape[0]-MALength):\n",
    "    mySlice = temp_X_Test[i:i+MALength]\n",
    "    MA[i,:] = np.average(mySlice,axis=0)\n",
    "\n",
    "temp = torch.from_numpy(MA)\n",
    "MAPredLoss = loss(temp,y_test[MALength:]).item()\n",
    "print(f'\\nThe loss on the TEST SET from using a moving average over the last 5 periods as a proxy for the factors at time t+1 gives us an average MSE of {MAPredLoss}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a37bbf0e-6006-4d22-b177-dc4a5309aa79",
   "metadata": {},
   "source": [
    "When we compare the trained model to both baseline cases we can see that we have managed to improve the Mean-Square Error. I will discuss the differences between the neural net from part 2 and the neural net from this part below. Before we do that, notice that the moving averge is set to 30 (representing an entire months worth of previous values). The moving average as a predictor actually does a fairly decent job compared to the model. I am suprised by the little improvement that the model has over this baseline method."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cc45d0f-ba10-4a1f-8329-79c267928b07",
   "metadata": {},
   "source": [
    "The main differences between this model and the previous model are the following:\n",
    "\n",
    "Learning Rate: The learning rate was increased from 0.01 to 0.03\n",
    "\n",
    "Number of Epochs: The number of epoch was finely tuned down to only 55, which is much less than before\n",
    "\n",
    "Sequence Length: The sequence length was adjusted to take in the previous 30 days worth of information (representing the previous month factor values)\n",
    "\n",
    "Hidden Dimension: The hidden dimension was increased from 2 to 5\n",
    "\n",
    "Dropout Rate: The dropout rate was decreased from 0.5 to 0.1 in this model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cee9930f-98bd-4d74-a92e-2c3ba98dec0d",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "Overall, I noticed that no matter what hyperparameters I changed I had trouble reducing the training and validation error below the current levels they are at. This was really interesting to observe. It seemed to be indicating a plateau in about the same region. I also noticed that it was much more difficult to overfit this data. I still managed to do it by cranking up the hyperparameter values to fairly large levels, but this model was MUCH more forgiving compared to model trained on monthly factors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "819b17cf-9420-471c-b189-a2dd9ddbb763",
   "metadata": {},
   "source": [
    "For the next part I am only going to repeat my strategy of buying 10,000 shares of the factor that shows the greatest t to t+1 difference and shorting 10,000 shares of the factor that shows the most negative difference from t to t+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1298,
   "id": "c9f4335a-a6a1-4e56-9b2d-11e66e81ea39",
   "metadata": {},
   "outputs": [],
   "source": [
    "predictedFactors = model(x_test)\n",
    "predictedDiff = ((predictedFactors - x_test[:,-1,:])).detach().numpy()\n",
    "actualDiff = ((y_test - x_test[:,-1,:])).detach().numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1299,
   "id": "2d0faf86-965e-4f81-a679-d22f4d4a4f94",
   "metadata": {},
   "outputs": [],
   "source": [
    "strategy = np.zeros(predictedDiff.shape)\n",
    "strategy[np.max(predictedDiff,axis=1).reshape(strategy.shape[0],1)==predictedDiff] = 1\n",
    "strategy[np.min(predictedDiff,axis=1).reshape(strategy.shape[0],1)==predictedDiff] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1300,
   "id": "098a5347-2025-48c1-a36f-aaf2637ddfb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "Diff1 = actualDiff*strategy\n",
    "Diff = np.sum(Diff1,axis=1).reshape(Diff1.shape[0],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1301,
   "id": "bf0dd7bb-c4fa-474a-a08a-d1203e9522a6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The total gain over the test set for buying 10,000 shares of the factor with the biggest expected difference and shorting 10,000 shares with the smallest expected difference is $12441649.98805616\n"
     ]
    }
   ],
   "source": [
    "totalNumberShares = 10000\n",
    "PNL = Diff*totalNumberShares\n",
    "print(f'The total gain over the test set for buying 10,000 shares of the factor with the biggest expected difference and shorting 10,000 shares with the smallest expected difference is ${np.sum(PNL)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee57fffd-4fdb-4e5e-ac02-83bdd18cea39",
   "metadata": {},
   "source": [
    "# Conclusion\n",
    "\n",
    "overall the daily optimized model is much more profitable than the monthly optimized model. What this seems to indicate is that we can predict the direction of the change in the Fama/French factors well enough on the daily level to capture more of the total price variation. I would be curious to see at what time aggregate level this strategy no longer works.\n",
    "\n",
    "Finally, it is fair to say that an appropriately optimized LSTM network does a decent job at improving the prediction of the t+1 Fama/French values as compared to using a simple moving average to predict the t+1 FF factor values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6954079-7b13-4547-96c0-4c4e44d96835",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0d3cbea-35a0-48f2-a78c-54fa5dd32a77",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}