{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extended Least Squares\n",
    "\n",
    "Example created by Wilson Rocha Lacerda Junior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we import the NARMAX model, the metric for model evaluation and the methods to generate sample data for tests. Also, we import pandas for specific usage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install sysidentpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sysidentpy.model_structure_selection import FROLS\n",
    "from sysidentpy.basis_function._basis_function import Polynomial\n",
    "from sysidentpy.metrics import root_relative_squared_error\n",
    "from sysidentpy.utils.generate_data import get_siso_data\n",
    "from sysidentpy.utils.display_results import results\n",
    "from sysidentpy.utils.plotting import plot_residues_correlation, plot_results\n",
    "from sysidentpy.residues.residues_correlation import compute_residues_autocorrelation, compute_cross_correlation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating 1 input 1 output sample data  \n",
    "\n",
    "The data is generated by simulating the following model:\n",
    "$y_k = 0.2y_{k-1} + 0.1y_{k-1}x_{k-1} + 0.9x_{k-2} + e_{k}$\n",
    "\n",
    "If *colored_noise* is set to True:\n",
    "\n",
    "$e_{k} = 0.8\\nu_{k-1} + \\nu_{k}$\n",
    "\n",
    "where $x$ is a uniformly distributed random variable and $\\nu$ is a gaussian distributed variable with $\\mu=0$ and $\\sigma$ is defined by the user.\n",
    "\n",
    "In the next example we will generate a data with 3000 samples with white noise and selecting 90% of the data to train the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train, x_valid, y_train, y_valid = get_siso_data(n=1000,\n",
    "                                                   colored_noise=True,\n",
    "                                                   sigma=0.2,\n",
    "                                                   train_percentage=90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the model\n",
    "\n",
    "First we will train a model without the Extended Least Squares Algorithm for comparison purpose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "basis_function = Polynomial(degree=2)\n",
    "\n",
    "model = FROLS(\n",
    "    order_selection=False,\n",
    "    n_terms=3,\n",
    "    extended_least_squares=False,\n",
    "    ylag=2, xlag=2,\n",
    "    info_criteria='aic',\n",
    "    estimator='least_squares',\n",
    "    basis_function=basis_function\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5649873632814074\n"
     ]
    }
   ],
   "source": [
    "model.fit(X=x_train, y=y_train)\n",
    "yhat = model.predict(X=x_valid, y=y_valid)\n",
    "rrse = root_relative_squared_error(y_valid, yhat)\n",
    "print(rrse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly we have something wrong with the obtained model. See the *basic_steps* notebook to compare the results obtained using the same data but without colored noise. But let take a look in whats is wrong."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Regressors   Parameters             ERR\n",
      "0    x1(k-2)   8.8975E-01  7.42588820E-01\n",
      "1     y(k-1)   3.0445E-01  8.42278843E-02\n",
      "2     y(k-2)  -4.7924E-02  2.09392505E-03\n"
     ]
    }
   ],
   "source": [
    "r = pd.DataFrame(\n",
    "    results(\n",
    "        model.final_model, model.theta, model.err,\n",
    "        model.n_terms, err_precision=8, dtype='sci'\n",
    "        ),\n",
    "    columns=['Regressors', 'Parameters', 'ERR'])\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Biased parameter estimation\n",
    "\n",
    "As we can observe above, the model structure is exact the same the one that generate the data. You can se that the ERR ordered the terms in the correct way. And this is an important note regarding the Error Reduction Ratio algorithm used here: __it is very robust to colored noise!!__ \n",
    "\n",
    "That is a great feature! However, although the structure is correct, the model *parameters* are not ok! Here we have a biased estimation! The real parameter for $y_{k-1}$ is $0.2$, not $0.3$.\n",
    "\n",
    "In this case, we are actually modeling using a NARX model, not a NARMAX. The MA part exists to allow a unbiased estimation of the parameters. To achieve a unbiased estimation of the parameters we have the Extend Least Squares algorithm. Remember, if the data have only white noise, NARX is fine. \n",
    "\n",
    "Before applying the Extended Least Squares Algorithm we will run several NARX models to check how different the estimated parameters are from the real ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\wilso\\AppData\\Local\\Temp/ipykernel_11404/3592060221.py:12: DeprecationWarning: setting an array element with a sequence. This was supported in some cases where the elements are arrays with a single element. For example `np.array([1, np.array([2])], dtype=int)`. In the future this will raise the same ValueError as `np.array([1, [2]], dtype=int)`.\n",
      "  parameters[:, i] = list(model.theta)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD7CAYAAACBiVhwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApYUlEQVR4nO3de3BU9fk/8Pc5e0mIxFvcFYqUaa0lLSrWG6U64avVAEqwBGYKtKJTL1AtKF/FC6SkVSlUsdGOZaY4VMfbt0WKBRkMWqmMGH/a8ocURKVKUkFIliDCJtns7jmf3x+75+zZ3YQsSc45n+S8X/+Q3cTdx5PsPvt8ns9FEUIIEBERWahuB0BERPJhciAiojxMDkRElIfJgYiI8jA5EBFRHiYHIiLKw+RARER5/G4H0F++/LINut73JRv33/+/AIDf/vZ3fX6s/nL//f+LQMCHRx55zO1QpFdWNhStrVG3w5Aer1NhZLtO/fn+pKoKzjjjlG6/P2iSg66LfkkOLS0t5uPJoqWlBcGgX6qYZMbrVBhep8LIdJ2cfH/isBIREeVhciAiojxMDkRElIfJgYiI8jA5EBFRHiYHIiLKw+QwgAghEN/1d4iYPPOuiWhwYnIYQPTIPnQ2vIDY//s/t0MhokHO1uQQjUYxZcoU7N+/P+v+F154ATfeeKN5e8+ePaiursbEiROxZMkSJJNJO8MasES8PfVv9IjLkRDRYGdbcvjggw8wa9YsNDY2Zt3/n//8B6tXr866b9GiRVi6dCm2bNkCIQTWrl1rV1gDm5pe0K5r7sZBRIOebclh7dq1qK2tRTgcNu+Lx+NYunQpFixYYN534MABxGIxXHTRRQCA6upq1NfX2xXWwCZ0tyMgIo+wbW+lZcuW5d33+OOPY/r06TjnnHPM+1paWhAKhczboVAIzc3NJ/18ZWVDexdojmDQn46jtF8erz8YMZ1WGkQHgEDAJ1V8suG1KQyvU2Fkuk5Ovj85tvHeO++8g4MHD+LBBx/Ee++9Z96v6zoURTFvCyGybheqtTXaL5tRxeOpfkckcrzPj9Vf4vEkgkE/vjqamqWUSOpSxSeTUKiU16YAvE6Fke069ef7k6oqJ/xQ7Vhy2LRpE/bu3YsbbrgB7e3tOHz4MO6++24sWrQIkUjE/LnDhw9nDUVRhmCvgYgc4lhyWL58ufn1e++9h6eeegpPPPEEAKCoqAg7duzAJZdcgg0bNqCiosKpsAYWMzmcfGVFRAOfEEAvBlZ6RYrzHFauXImamhpEo1GMGTMGc+bMcTskOelsSBN52SefH8UpQwKOPJftyWHr1q15940bNw7jxo0zb5eXl2PdunV2hzLwGZWDUx8diEg6bR0JR56HK6QHEKGnFwcyORCRzZgcBhL2HIg8y+njSpkcBhL2HIg8K550drYik8NAwqmsRJ4VTzr74ZDJYQARbEgTeVaSyYG6ZTak+Wsj8hpdsOdA3TE23mPlQOQ5DucGJocBhT0HIs8SluzgRBXB5DCQMDkQeZY1HyQc6D8wOQwgZkPa6fqSiFxnrRY6E/Z/UGRyGEjM5MD1DkReY/1MGGdyoCysHIg8S2RVDhxWIiuz58DkQOQ11lc9KwfKIgQrByKv4rASdY89ByLP4rASdY89ByLPYuVA3TN2ZWVyIPIcTmWlbhmH/Qg2pIk8h5UDdY/DSkSeZe05dMSZHMiKyYHIs6wv+7aY/edI25ocotEopkyZgv379wMA/vKXv2DKlCmoqqrCgw8+iHg8DgDYs2cPqqurMXHiRCxZsgTJZNLOsAYuc5YSkwOR11h7DvH4AJ6t9MEHH2DWrFlobGwEAOzbtw9r1qzBn//8Z2zcuBG6ruOll14CACxatAhLly7Fli1bIITA2rVr7QprQBOcykrkWdZhJc2BI4NtSw5r165FbW0twuEwACAYDKK2thZDhw6Foij49re/jS+++AIHDhxALBbDRRddBACorq5GfX29XWH16MvOr7D36GfQZXwD5rASkWdZX/ZJzf73AL9dD7xs2bKs2yNGjMCIESMAAEeOHMGLL76I5cuXo6WlBaFQyPy5UCiE5ubmk36+srKhfQs47UjsCABgyKkqTi0u7ZfH7KtgMPVr8qtAHIDfpyAUkiM2GfHaFIbXqTCyXKcvjsbMr/0Bn+1x2ZYcutPc3Ixbb70V06dPx7hx47Bjxw4olpPNhBBZtwvV2hqFrvc9mxrZ+YuWI+gc0ueH6xfxeBLBoB/JdI8mmdAQiRx3OSo5hUKlvDYF4HUqjEzX6cuj7ebX0fZ4n+NSVeWEH6odna306aefYubMmZg2bRruvPNOAMCwYcMQiUTMnzl8+LA5FOUGYw1BUpevKW72HCDhkBcR2Sqr56AN4J5Drmg0iltuuQV33XUXfvazn5n3jxgxAkVFRdixYwcAYMOGDaioqHAqrG4lhYSnrrHnQORZ1pe91g+jJD1xbFhp3bp1OHz4MJ555hk888wzAICrr74ad911F1auXImamhpEo1GMGTMGc+bMcSqsbslYOTA5EHmXWTkoQNKBysH25LB161YAwM0334ybb765y58pLy/HunXr7A7lpCRlPK85HRO3zyDyHj2TG/qlv9oTrpDOoSDVDGflQEQyMSoHRVHgQG5gcuiOjD0HweRA5FmWUaWs5rRdmBy6oclYOQiukCbyKmvPQWdycI8u47i+4HkORF5l5gZFgQO7ZzA5dEfO7TOMpMDkQOQ1RrXAYSWXSZkcwMqByKvMlz2HldwlZXIw/iBkjI2IbCUslQOHlVwkW3LIKiNZORB5jrXn4MRaJyaHbsiWHIjI2/SsyoHJwTXyJYfMH4OQLjYislt2z8H+52Ny6IZ0ycH6x8BhJSLPyfQcFM5WcpMTswF6TebYiMgWZuGgcFjJVbp022eIbr4mIi/IXudg//MxOXRDyhXSAKD4OJWVyIOEpXTgOgcX6U5MJD4Zxh+DqnJYiciDstY5MDm4R5f1KE7Vx+RA5EGZdQ4cVnJJ6qprsg7dqBxWIvIi3borKxvSzjO3tpPtDdg86EMFG9JE3pPJDZzK6jjr2ga5Kwe3gyAip2VOguMiOMdZs7ETmflkmHuppCsH2eIjIntZX/IDviEdjUYxZcoU7N+/HwDQ0NCAqqoqVFZWoq6uzvy5PXv2oLq6GhMnTsSSJUuQTLpzCpuWVTnIts4hTfWlv2ByIPIS6xnSYiD3HD744APMmjULjY2NAIBYLIbFixdj1apV2Lx5M3bt2oVt27YBABYtWoSlS5diy5YtEEJg7dq1doV1QtaFb9J9MjfCMZKDbPERka30wTKVde3ataitrUU4HAYA7Ny5E6NGjcLIkSPh9/tRVVWF+vp6HDhwALFYDBdddBEAoLq6GvX19XaFdULWnoN8U1nTfxhq+lfG5EDkKU5PZfXb9cDLli3Lut3S0oJQKGTeDofDaG5uzrs/FAqhubn5pJ+vrGxo74NNOxbLfB0s8iEUKu3zY/aHYNAPPZGqavyBIOIAzjrrFKj+oLuBSUqW35vseJ0KI8t1KilJvd59PhVCCNvjsi055NJ1HYqimLeFEOmDsru+/2S1tkb7PPf3q85j5tftHZ2IRI736fH6SzyehD9dOSTTBc3hyDEo/iIXo5JTKFQqze9NZrxOhZHpOkWjnQBSaxx0gT7HparKCT9UOzZbadiwYYhEIubtSCSCcDicd//hw4fNoSinDYyprBxWIvKi3D6D3X1Rx5LD2LFjsW/fPjQ1NUHTNGzatAkVFRUYMWIEioqKsGPHDgDAhg0bUFFR4VRYWXSJp7IaFIUNaSIvMwZW7H4LcGxYqaioCCtWrMD8+fPR2dmJCRMmYNKkSQCAlStXoqamBtFoFGPGjMGcOXOcCiuLLvNUVuvGe6k7XAuFiJynC6SmKpm3BVSc/BB8oWxPDlu3bjW/Hj9+PDZu3Jj3M+Xl5Vi3bp3dofTIOkNJ1soBavpXJmt8RGQLIQQUZPLDoBlWGgiE1D2HnMqByYHIUzIv+VR6sHsdHJODhbXnIN0Z0oZ0z0G6jQGJyFZmpaDk3LYJk4NF1iI42d58za2V2JAm8iLzILj0v3afR8bkYGHtOci2QlrkDiuxIU3kKblrwITN7wFMDhbWakHahjSnshJ5Uu5L3u63ACYHC6l7DrlTWZkciDzFurcSYP9pcEwOFlL3HNIUbtlN5EmZ0Qwl57Y9mBwshNTJwfjYYAwryRYfEdkpZw0cp7I6KbNfuuLIfum9wmElIk+yHhNqvW0XJgcLo1pQFEW62Uo87IfI2wSQt32GnZgcLIyEkDppSbLkYE5lZXIg8qKcNXCcreQks3KAImFySDEa0kK2yoaI7DVYt+weCMyeg6JIuM7BXCKdvilbfERkJ12kqwaHtuwuKDnMnz8fDQ0N9kYiASFz5WDOYlOzbxORR2S/6KXoOVx77bVYtWoVJk6ciDVr1uDo0aO2BuUWHZnKQbqGtMFc5yBpfERki9R5DgoUmXZlnTp1Kl544QWsWrUKra2tmDFjBhYtWoSdO3faG53DdD11wI+C1NnWMhFsSBN5W972GRJUDgCg6zqamprQ2NgITdNQVlaGX/3qV/j9739vZ3yOyq4cJHvz5a6sRJ5mHvYj0zGhdXV1WL9+PUaOHInZs2fjySefRCAQQHt7O6666iosWLDA3igdIvf2GawciLwsdxjJ7r2VCkoOR44cwdNPP43y8vKs+0tKSvD444/bEpgbrCukpT1Mx2xISxofEdkkZyqrDFt2a5qWlxiMauHKK6/s/6hcIiwrpKU7JjR3hbRsw15EZCunF8GdsHKora1Fc3MzduzYgSNHjpj3J5NJfP75571+0g0bNmD16tUAgIqKCtx///1oaGjA8uXL0dnZicmTJ2PhwoW9fvzeyqyQlnedA3sORN6UnqxkNh3snsp6wuQwY8YM7N27Fx9//DEmTpxo3u/z+XDRRRf16gk7OjqwbNky1NfX49RTT8WsWbOwdetWPPTQQ3j++ecxfPhwzJ07F9u2bcOECRN69Ry9pcm8t5IhPawkX/IiIjvlvubtHtw4YXK44IILcMEFF+CKK67A2Wef3S9PqGkadF1HR0cHSkpKkEwmMXToUIwaNQojR44EAFRVVaG+vt7x5CCydmWVLTnkNKQ5rETkKbnDSq5WDnfddReefPJJ3HrrrV1+/9VXXz3pJxw6dCjuuusuTJ48GUOGDMFll12GlpYWhEIh82fC4TCam5tP6nHLyoaedCy5So4EAAA+nwodAqFQaZ8fsz8Eg34kY6mvTz39FMQAnH5aMYZIEp9sZPm9yY7XqTCyXKdg0A9FUeD3p0YPTjttiK2xnTA53HbbbQCAX/7yl/32hB999BH++te/4h//+AdKS0tx7733orGxMfvg7JyDtAvR2hrt89SuY8c7Us+vp9Z1RCLH+/R4/SUeT5ozB45F4wCAo1+2ITpEjvhkEgqVSvN7kxmvU2Fkuk6xzgQggKSmIwDgyJftfYpNVZUTfqg+4Wyl888/HwBw+eWXY/jw4bj88svR3t6Of/7zn/jOd77Tq4C2b9+O8ePHo6ysDMFgENXV1XjvvfcQiUTMn4lEIgiHw716/L7IDCvJu85BUdiQJvIikT4KLjNbSYKprEuXLsXTTz+NTz/9FDU1Ndi/fz8WL17cqycsLy9HQ0MD2tvbIYTA1q1bMXbsWOzbtw9NTU3QNA2bNm1CRUVFrx6/L6wNaQEhV9OXh/0QeVpeQ1qGFdK7du3CunXrsHr1akybNg333HMPqqure/WEV155JT788ENUV1cjEAjgggsuwPz583HFFVdg/vz56OzsxIQJEzBp0qRePX5fWKeyAqnqwWd8UncdG9JEXmYmg/SQu90fXgtKDkIIqKqKd955B/PmzQMAxGKxXj/p7bffjttvvz3rvvHjx2Pjxo29fsz+kFkVnU4OEJAsNXCFNJHHOTVbqaBhpa9//eu47bbbsH//flx++eW45557MHr0aFsDc4Oet7GVhG/AHFYi8iTj/Slz297nK6hyWL58Od544w1ccsklCAQCuPTSS/GjH/3I3shckGpCK4BlWEkaRrNcNfI5kwORF0nVkC4pKcGll16KY8eOYffu3bjwwgvx2Wef2RqYG3ToltQgWXIwcLYSkSfpQlg/u8rRkH7yySfxpz/9CWVlZeZ9iqLgzTfftC0wN+hCNzYvSd+W6Q1YpGJzqBlFRHLJfclLsWX3hg0b8Prrr/fbFhqy0i3rHADItb+SAADVuZM+iEguRk/UuGnz0xU0rDR8+PBBnxiAVANaAcyrL9ewkgBUxZIcZIqNiOyWKhQy6UGKqazjx4/Ho48+ih/+8IcoLi427x8zZoxtgbnBaEgrMjakAaRyORvSRJ5mfniVIDmsX78eAFBfX2/eNzh7DrlTxSR7A1Y5rETkVSJn2NvVLbsNW7dutTcKSehINaSlrBzyN1ZxMxoicpgQmc+GgCSL4Nra2vDQQw/hpptuwtGjR7F06VK0tbXZGpgbdKPnkCbdIjhFgcJhJSJPcnpvpYKSwyOPPILS0lK0traiqKgI0WgUS5cutTcyF5g9h3R6luscaQFF4bASkVflbK0kxyK4PXv2YOHChfD7/RgyZAhWrlyJPXv22BqYG0ROz0FI9OlcGAfIMjkQeVLeOgcZkoOqZv+Ypml59w0GRo9B2hXSCqeyEnmVgPHh1ZjKau/zFdSQvuyyy/DYY48hFovh7bffxgsvvIBx48bZG5kLNKGnh5QkbEhDpHZkNXdlZeVA5CXmnBSZhpXuvfdelJSUoLS0FE888QTKy8tx33332RqYG3ShIdVzMG5L9AYskE4M6U8NEg15EZH98oeV7H2+HiuHN954A2vWrMHHH3+M4uJijB49GhdffDGKiorsjcwFmjlbSdbKgT0HIq8SRumQ5uoiuNdeew11dXVYsGABysvLoSgK/v3vf2PZsmXo7OxEZWWlrcE5TdM1c6YSIFtyQHpYiT0HIi8yZysZt92sHJ577jk8++yz+NrXvmbed+6552Ls2LFYvHjx4EsOQktvnmEM3Uj2BmzphxCRt5izKWXoObS1tWUlBsM3vvENdHZ22haUW3IrBbl6DoKVA5GHSTWV1efr/gTlwXiegCZSw0pyLoJDOjbOViLyoswOOs5MZXVlscLWrVtRXV2NyZMn45FHHgEANDQ0oKqqCpWVlairq3MjrHRDOjNsI9f2GbmVA5MDkZfkzlB0dcvujz/+GBdffHHe/UIIxOPxXj3h559/jtraWrz88ssoKyvDTTfdhG3btqG2thbPP/88hg8fjrlz52Lbtm2YMGFCr56jtzRdAyDpIjiukCbytvQpoZktu+19uhMmhzfeeKPfn/CNN97Addddh2HDhgEA6urq0NTUhFGjRmHkyJEAgKqqKtTX1zufHISWvQhOtrUEWQ1piRIXEdnOelKlACDcPCZ0xIgR/f6ETU1NCAQCmDdvHg4ePIj/+Z//wXnnnYdQKGT+TDgcRnNz80k9blnZ0L4HpwI+VUUgkOq1DC0NIhQq7fvj9lEw6EccAv5AAGeFTkUUwCklQZwuQWwykuF3NhDwOhVGluvk86lQVQWBgA9xAENK7H1/Kmj7jP6kaRr+9a9/4fnnn0dJSQl+/vOfo7i4OGt9gRAi63YhWlujfT5wO5FMQtcFksnUp/KjX7UjUny8T4/ZH+LxJASApCZwuDUKAIhGY0hE3I9NNqFQKSK8Lj3idSqMTNcpkdShC4FEQoMCIBrt7FNsqqqc8EO148nhrLPOwvjx43HmmWcCAK655hrU19dnzYyKRCIIh8NOh5ZaBGdpScvVcxA5w0qSDXkRkb2EMN+dFEWxfQsdx2crXXXVVdi+fTuOHTsGTdPw9ttvY9KkSdi3bx+ampqgaRo2bdqEiooKp0NL9xwApw7wPmnceI/Is6wveUUBdBmOCe1PY8eOxa233orZs2cjkUjgiiuuwKxZs/DNb34T8+fPR2dnJyZMmIBJkyY5HVrmsJ/0benWOYCzlYi8ylopqKri7lRWu8yYMQMzZszIum/8+PHYuHGjG+GYjO0zDNKtc1C5zoHIq3Irh0G5CE5WqfMcYDbDdZmmiwogNdzF7TOIvEiITMdRURQ5ToLzAiEE9JwV0lLtrZSuHDLrMGSKjYjsJpDJDqoiyTGhXqAJLf1VZr6SVLOVAJh/GU7UlEQkFetLXlUUDis5xTw/Wsm/TwrG9hkAACYHIq+xNqAVxf6GNJNDmlE5KJZxfamSg7HxHpCuHGSKjYjsljeVlZWDMzRz0rBimRAk2adzRcn+l4g8w3pIqMrKwTlm5WDZL12qdQ5CZM5yUBTJptkSkd2EuUuC0XZkcnBEUk8CQPZsJZmmsgKWYSWVPQcij8keVlI4rOSUhJ4AgJwNAOVJDqlpbNaGtDyxEZH9rJWCqti/ZTeTQ1oiXTmo6cpBVVTJ1jkgUzmoqv0bqxCRdLIXwdn7XEwOaUZyMMb1U8lBojdgy1RWRfUBunbinyeiQUXPXecw2HZllVXSGFYyKgcociUH61RW1Q/B5EDkKdnrHLi3kmPiWnbPQVVUCRvSxtp5H5CudIjIGyyTlVLDSuw5OCN3tpIiW8/B+pfBYSUiz7G+HzmxZTeTQ5rZkE6/Aftk6zkg0w9hz4HIe1KVgnWdg73Px+SQlsirHGTrOSCrcmDPgchbsoaVwC27HZO7zkGFKtU6h+yGtA8QTA5EXqLnrnNg5eCMRO5sJel6DshODqwciDzD6C+Y6xxUVg6OSZo9B8s6B6lmK1nGG5kciDzFTASWw35YOTgkoSez9lVSZew5qFznQORFuRsiDPpjQn/729/igQceAAA0NDSgqqoKlZWVqKurczyWhJZAQPWbt2VbIS04rETkWZlhJQ/syvruu+/ilVdeAQDEYjEsXrwYq1atwubNm7Fr1y5s27bN0XgSehIBNWDelm6dg2XjPUX1ARoXwRF5Rf6w0iA9JvTo0aOoq6vDvHnzAAA7d+7EqFGjMHLkSPj9flRVVaG+vt7RmJJ6An5L5eBTZJutlFnnwMqByFtycoMjx4T6e/6R/rd06VIsXLgQBw8eBAC0tLQgFAqZ3w+Hw2hubj6pxywrG9qnmNRPFRQHixAMpi5J0O+HL6giFCrt0+P2h2DQjw4hUHJKMc4MlaK5pBidR3UpYpMRr0theJ0KI8N1irbHAQA+n4pg0I+ioB+6ELbG5nhyePnllzF8+HCMHz8e69evBwDoup5zjoLIul2I1tZon/Yaiba3QxUq4vHUcI2mCcQ644hEjvf6MfuLEVN7RwJa5Dg64wJaIiFFbLIJhUp5XQrA61QYWa7T8XRy0DQd8XgSyaSGhKb3KTZVVU74odrx5LB582ZEIhHccMMN+Oqrr9De3o4DBw7A5/OZPxOJRBAOhx2NK9VzyG5Iy3WGdGYRnOLzs+dA5CHmsFLWVNZBNqz0zDPPmF+vX78e77//Pn7961+jsrISTU1NOOecc7Bp0yZMnz7d0bgSehJ+S0NattlKADKzlXx+CO7KSuQZmckxme197P7s6krPIVdRURFWrFiB+fPno7OzExMmTMCkSZMcjSGpJxBUg+ZtudY55Hxs8AWA9BbjRDT4ddWQtnvLbleTQ3V1NaqrqwEA48ePx8aNG12LJaElUOIvMW+rUKHJtn+RuSsrh5WIvMRMBFwh7byueg7SrHMwK0pL5SB0CJ4jTeQJue9FTkxlZXJI67LnIM3eSunVkZaeAwBA59ASkRfkbbynZJ8pbQcmh7SknkDQZ60cZOo5pFlnKwEcWiLyCD27H51eIc3KwRHxnMpBkWm2Uu48Nl8qTsGmNJEnGD0Hc28lbtntnKSev/GeXOsckN2QBlg5EHlE7nsRG9IOEULkN6Qh4bCSmtNzYHIg8oTc/gIb0g5JpqesBqRdBJcz4GgMK7EhTeQJWnpmonmGtJI/g6m/MTkgNaQEQN7zHIy/AZUNaSIvSmrGjEXjGONBumW3bOLpN9n8qayy9Bxyp7IaDWkmByIv0LR05ZC+PagP+5FJV5WDIuVUVmMem1E5cFiJyAuSxmwlc1hJ4ToHJyTSm9gFfJnKwSfTsJK5txKHlYi8yKwczGEl9hwcYSaHrMpBouRg9qNzhpXYkCbyBE3LXiHNnoNDEuk32ayeA2Rc52DMVmLlQOQlyZyN9ziV1SFdz1ZSpNtbKTOslE5iTA5EnpDMGVZSFNi+ZTeTA7oeVlIVFZosw0oGY1gpHSe3zyDyBg4ruSSTHKwNaR90XZLzHPL2VuKwEpGXJPXcykGBsHmqPZMDUgf9ANmVg1/1Iyk0qfoOSu6wEhvSRJ5gVg5ZK6TtfU4mB2QqB2tD2p9OFEmZToPLqRy4CI7IG3IXwamKAsGeg/2MhnTQZ00OvvT3JHgDFtkNaXBXViJPySyCM7bsZuXgiEzlkD2sBMiRHMy/AWNYSVFSCYINaSJPyCyCS90etIf9PPXUU7j++utx/fXX49FHHwUANDQ0oKqqCpWVlairq3M0nkQXU1n9ikSVA3Ia0gDg80PI0jAnIlsZG+8ZlMF4nkNDQwO2b9+OV155BX/729+we/dubNq0CYsXL8aqVauwefNm7Nq1C9u2bXMspoSehAIFvnRCAKyVg0RvwIpq+ZKVA5FXJHUdfl/mw6GCQXgSXCgUwgMPPIBgMIhAIIBzzz0XjY2NGDVqFEaOHAm/34+qqirU19c7FlMifQqcYvlknmlIS1A55GzZDSDVlJaiqiEiu2magM/y+lfV1HuVnQnC3/OP9K/zzjvP/LqxsRGvvfYafvrTnyIUCpn3h8NhNDc3n9TjlpUN7XVMgf+qCPqDCIVKEQymLknZ6aUAgFNPK0LojNJeP3Z/8PsVxAGcccYpKA6lYukIBFEUUBEKuRubjHhNCsPrVBgZrlMg6EcwoALp96fS0iIAwJlnDkXAb89nfMeTg2Hv3r2YO3cu7rvvPvh8PjQ2NprfE0JkfYovRGtrtNfLyY+1tcMHHyKR44jHU5/G247HAQAtrcdwSvJ4rx63vySTqWbU0a9i8BWlYtGgItbWjkjE3dhkEwqV8poUgNepMLJcp+PRGFRVMd+f4rHUv4eav0JxsHdv46qqnPBDtSsN6R07duDmm2/GPffcg2nTpmHYsGGIRCLm9yORCMLhsGPxGMNKVjLNVsrdWyn1pR+QqR9CRLZJJAUCvszr35fuP+Q2qvuT48nh4MGDuPPOO7Fy5Upcf/31AICxY8di3759aGpqgqZp2LRpEyoqKhyLKaEns85yACTtOSjZPQfurUTkDUlNh9+SHIyvjSmudnB8WGnNmjXo7OzEihUrzPtmzpyJFStWYP78+ejs7MSECRMwadIkx2JKdlE5BGSsHNTMbKpUQ5qVA5EX5CUH1f7KwfHkUFNTg5qami6/t3HjRoejSUnoyaytMwDJprIaMxLU3GElGRIXEdktoekI+C2zKdOJwtiQzw5cIY2uew7GDq1xLe5GSF1SrDFyWInIM5JJHb7B3nOQUUJPZm3XDQBD/MUAgJjW6UZI2XL3VgJSQ0yaBFUNEdkuqWU3pJ3oOTA5wEgO2ZVDkS81jziWjLkRUo78noPiC3BYicgjEnkNaVYOjkhqibyeQ0D1Q1VUOSqHrhrSqp9bdhN5hKZlb59hDDElWTnYK6EnEfRlVw6KoqDYV4ROGZKDOapkna3kY+VA5BEJTWSthDZmK3FYyWYJPb9yAIBifzFiSQmSQ5eL4ALceI/II5LJrtc5JDisZK+ueg4AUOwrkmNYSXS1zsHHLbuJPMKNRXCeTw5CCCS7Sw7+Ikka0mlZySHAk+CIPCKp6V1vn2HjcXCeTw6a0CAg8qayAqkZS1IMK3UxlTV1ngOTA5EXJDQd/q4WwbFysE9Xp8AZSvxD0KF1OB1SFwQAJXunWtUHCA1C2PfHQURySCZFN9tnMDnYJq4Z50fnVw4lgRK0J9xPDkKI7CNCgdSwEsD9lYgGOV0X0EXOIji/0XPgsJJtjO0xinzBvO+dEihBW6IdugyfznOSg+JL9x84tEQ0qCXS1YHf38VspSQrB9sY6xiK/EV53zvFPwQCwv2+gxDIO/ooPQwmuNaBaFAzho6sw0pFgdSHw86EfSMHTA4nqBxKAiUAgPZku6Mx5RPZ+yoBmWElVg5Eg5pxEqR1hbTfp8CnKkwOdjIrB18XlUM6ObQlXE4OQs8fVlI5rETkBcbQkXWFtKIoKAr4EIszOdjmhJWDP105uN2UFgJKd5UDh5WIBrVYujowhpIMRUEfKwc7nbhyGAIAaEu0ORpTHqF3O6wkkvKcN0FE/a+zu+QQ8KGTlYN9YmZy6Gq20ikAgLak+5VD3rBSMJW4RNz9qbZEZJ94nJWDK6LxVFVQ4h+S9z3jvnbXew75DWmlKDXkJeJuN8uJyE6diVTPoSiYnRyK2XOwVzTRhlMCJfCpvrzv+VQfin1FOC7BsJKi5lYOqeQAJgeiQS2WSPUVcyuH4qAPsU77eo5SJYdXX30V1113HSorK/Hiiy868pzReBRDA0O7/X7ZkDPR2tHqSCzdEbpmrmswGMlBdDI5EA1m0fbUFj9DS7J3cThtaBGORu1bg5W/oZBLmpubUVdXh/Xr1yMYDGLmzJkYN24cvvWtb9n6vJGOVpxZfHq33w+XhHDg+Be2xnAiQtcBXcusiDYEhwCKD6LjK3cCIyJHHGuPQ1UUDB2SnRzOLC3CsfYEEkk9a5prf5EmOTQ0NOD73/8+Tj/9dADAxIkTUV9fj1/84hcF/feqmreGuEfH4seREHF856zzzP8+HA5nPd53zvwW9kf3IxI7jLNLQif9HH2lH4sgPPxrCJxxdvb/o+pH8OtjgPYjvfp/H8x4PQrD61QYt69TtCOJ0V8/HX6fmvX+9PXhpQifMQQtRzswMtz96Ed3evr/UoQQ9u3cdBL++Mc/or29HQsXLgQAvPzyy9i5cycefvhhlyMjIvIeaXoOuq5nbUkthMjeopqIiBwjTXIYNmwYIpGIeTsSiZglFBEROUua5PCDH/wA7777Lo4cOYKOjg68/vrrqKiocDssIiJPkqYhffbZZ2PhwoWYM2cOEokEZsyYgQsvvNDtsIiIPEmahjQREclDmmElIiKSB5MDERHlYXIgIqI8TA5ERJSHyUEyPW0++Pe//x033HADpk6dijvuuANffeXNvZUK3aTxrbfewtVXX+1gZHLp6Tp99tlnuPHGGzF16lTccsst/Hvq5jrt3r0b06dPx9SpUzF37lwcO3bMhSgdJkgahw4dEldddZX48ssvRVtbm6iqqhJ79+41v3/8+HFxxRVXiEOHDgkhhHjiiSfEww8/7Fa4runpOhkikYiYNGmSuOqqq1yI0n09XSdd10VlZaXYtm2bEEKIxx57TDz66KNuheuaQv6eZs2aJd566y0hhBDLly8Xv/vd79wI1VGsHCRi3XywpKTE3HzQkEgkUFtbi7PPPhsAMHr0aBw8eNCtcF3T03Uy1NTUFLxx42DU03XavXs3SkpKzMWm8+bNw09+8hO3wnVNIX9Puq6jrS11rktHRweKi4vdCNVRTA4SaWlpQSiU2fk1HA6jubnZvH3GGWfg2muvBQDEYjGsXr0a11xzjeNxuq2n6wQAzz33HL773e9i7NixTocnjZ6u03//+1+cddZZWLx4MaZNm4ba2lqUlJS4EaqrCvl7euCBB1BTU4Mrr7wSDQ0NmDlzptNhOo7JQSKFbj54/Phx3H777SgvL8e0adOcDFEKPV2nTz75BK+//jruuOMON8KTRk/XKZlM4v3338esWbPwyiuvYOTIkVixYoUbobqqp+sUi8WwZMkSPPvss9i+fTtmz56N+++/341QHcXkIJFCNh9saWnB7NmzMXr0aCxbtszpEKXQ03Wqr69HJBLB9OnTcfvtt5vXzGt6uk6hUAijRo3CBRdcAACYMmUKdu7c6XicbuvpOn3yyScoKioyt/P58Y9/jPfff9/xOJ3G5CCRnjYf1DQN8+bNw+TJk7FkyRLPbmne03VasGABtmzZgg0bNmD16tUIh8N46aWXXIzYHT1dp+9973s4cuQIPvroIwDA1q1bMWbMGLfCdU1P12nUqFE4dOgQPvvsMwDAm2++aSbUwUyajfeo+80Hb7vtNixYsACHDh3Chx9+CE3TsGXLFgDA+eef77kKoqfr5IUXbiEKuU5/+MMfUFNTg46ODgwbNgyPPvqo22E7rpDrtHz5ctx9990QQqCsrAy/+c1v3A7bdtx4j4iI8nBYiYiI8jA5EBFRHiYHIiLKw+RARER5mByIiCgPkwMREeVhciAiojxMDkRElOf/A41pv39A+9PEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "parameters = np.zeros([3, 50])\n",
    "\n",
    "for i in range(50):\n",
    "    x_train, x_valid, y_train, y_valid = get_siso_data(n=3000,\n",
    "                                                       colored_noise=True,\n",
    "                                                       train_percentage=90)\n",
    "    \n",
    "    model.fit(X=x_train, y=y_train)\n",
    "    parameters[:, i] = list(model.theta)\n",
    "\n",
    "sns.set()\n",
    "pal = sns.cubehelix_palette(3, rot=-.5, dark=.3)\n",
    "\n",
    "ax = sns.kdeplot(parameters.T[:, 0])\n",
    "ax = sns.kdeplot(parameters.T[:, 1])\n",
    "ax = sns.kdeplot(parameters.T[:, 2])\n",
    "# plotting a vertical line where the real values must lie\n",
    "ax = plt.axvline(x=0.1, c='k')\n",
    "ax = plt.axvline(x=0.2, c='k')\n",
    "ax = plt.axvline(x=0.9, c='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the Extended Least Squares algorithm\n",
    "\n",
    "As shown in figure above, we have a problem to estimate the parameter for $y_{k-1}$. Now we will use the Extended Least Squares Algorithm.\n",
    "\n",
    "In SysIdentPy, just set *extended_least_squares* to *True* and the algorithm will be applied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\wilso\\AppData\\Local\\Temp/ipykernel_11404/3433714249.py:20: DeprecationWarning: setting an array element with a sequence. This was supported in some cases where the elements are arrays with a single element. For example `np.array([1, np.array([2])], dtype=int)`. In the future this will raise the same ValueError as `np.array([1, [2]], dtype=int)`.\n",
      "  parameters[:, i] = list(model.theta)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD9CAYAAABX0LttAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhhklEQVR4nO3dfXBU1f0/8Pe9u3kgJSqNu4amKd+WUmIRAaUq1QlDWx4sBAWZloeBdlolVoVKLVRDJD5hGEyLMmpbHTqOD52RUjHKlCCVytTGGdr8gYVS5CeGaoCwJJDn3c3ee35/7N6bfZBkIbn3ns19vzod3Zuw+/G47mc/53POuYoQQoCIiCiO6nQAREQkHyYHIiJKweRAREQpmByIiCgFkwMREaVgciAiohReK5/82Wefxe7duwEA06dPx7p16/DQQw+hoaEBI0aMAADcd999mDlzJo4cOYL169ejq6sLU6dOxaOPPgqv19LwiIjoAiz79K2vr8f777+PnTt3QlEU3Hnnndi7dy8OHTqEV199FX6/P+H3165diyeeeAKTJ09GRUUFtm/fjqVLl1oVHhER9cOyaSWfz4cHH3wQ2dnZyMrKwtixY3Hy5EmcPHkSFRUVKCsrw9atW6HrOpqamhAMBjF58mQAwMKFC1FXV2dVaERENADLKodx48aZf9/Y2Ijdu3fjtddew4EDB1BVVYX8/HyUl5djx44dGDduHHw+n/n7Pp8Pzc3NVoXWr5/97GcAgN/+9reOvP7nkTEmInKGXZ8Hlk/qHzt2DOXl5Vi3bh2+9rWv4bnnnjN/tnz5crz55psYO3YsFEUxrwshEh6no6WlE7o++JNAmppOAQACgY5BP9dQaWo6hexsr1Qxyczny+dYpYHjlB7ZxmmoPqNUVUFBwcgL/3xQzz6AhoYG/PjHP8YDDzyABQsW4OjRo9izZ4/5cyEEvF4vCgsLEQgEzOtnz55N6UkQEZF9LEsOp06dwr333ouamhrMnTsXQDQZPPnkk2hra0Nvby9ef/11zJw5E0VFRcjJyUFDQwMAoLa2FqWlpVaFRkREA7BsWmnbtm0IhULYtGmTeW3x4sVYuXIllixZgkgkglmzZmHevHkAgJqaGlRWVqKzsxMTJkzAihUrrAqNiIgGYFlyqKysRGVl5ef+bNmyZSnXSkpKsGPHDqvCISKii8Ad0kRElILJgYiIUjA5EBFlECGA851hWH0TTyYHIqIMEjjfg+bWbhz8fy2Wvg6TAxFRBtFim327gr2Wvg6TAxFRBjEOj9CG4ESI/jA5EBFRCiYHIqIMxIY0ERHZjsmBiIhSMDkQEVEKJgciIkrB5EBElIGsbUczORARZSaLswOTAxERpWByICLKIIpNr8PkQEREKZgciIgoBZMDEVEG4molIiKyHZMDERGlYHIgIspEPJWViIhMij2LWZkciIgyEBvSRERkOyYHIqIMZPXkEpMDERGlYHIgIsoksVVK7DkQEVEKi1eyMjkQEWUiwX0ORESUjJUDERGlYM+BiIhScFqJiIhSZPS00rPPPou5c+di7ty52Lx5MwCgvr4eZWVlmDVrFrZs2WL+7pEjR7Bw4ULMnj0b69evRyQSsTI0IqKMJiyeWLIsOdTX1+P999/Hzp078eabb+Lw4cPYtWsXKioq8Pzzz+Mvf/kLDh06hP379wMA1q5diw0bNmDPnj0QQmD79u1WhUZElPEytnLw+Xx48MEHkZ2djaysLIwdOxaNjY0YM2YMiouL4fV6UVZWhrq6OjQ1NSEYDGLy5MkAgIULF6Kurs6q0IiIMl7G9hzGjRtnftg3NjZi9+7dUBQFPp/P/B2/34/m5macOXMm4brP50Nzc7NVoRERZTyrKwevtU8PHDt2DOXl5Vi3bh08Hg8aGxvNnwkhoCgKdF2HEndGuXH9YhQUjBySeLOzo0Pi8+UPyfMNBRljkh3HKj0cp/TINE6qGv1On5eXbWlcliaHhoYGrF69GhUVFZg7dy4OHDiAQCBg/jwQCMDv96OwsDDh+tmzZ+H3+y/qtVpaOqHrg0+l4XAkFlvHoJ9rqITDEWRne6WKSWY+Xz7HKg0cp/TINk6argMAOjtDg4pLVZV+v1RbNq106tQp3HvvvaipqcHcuXMBAJMmTcInn3yCEydOQNM07Nq1C6WlpSgqKkJOTg4aGhoAALW1tSgtLbUqNCKijGf1JjjLKodt27YhFAph06ZN5rXFixdj06ZNWLVqFUKhEKZPn445c+YAAGpqalBZWYnOzk5MmDABK1assCo0IqKMZ3VD2rLkUFlZicrKys/92VtvvZVyraSkBDt27LAqHCKiYWUIZtH7xR3SREQZKUOXshIRkQViOSFjN8EREdHQM3ICkwMREaXI2B3SRERkHVYORERkMk5jzdhTWYmIyAJsSBMR0YXo7DkQEVEKVg5ERGTgUlYiIkoVSwo6G9JERJSClQMRERn6ppVYORARURKeykpERHFE0l+tweRARJRJjIY0KwciIjKw50BERBfGyoGIiJJxWomIiEx97WhOKxERkcGexUpMDkREmYinshIRUZzYzX5YORARkcGmWSUmByKiTMR9DkRE1Ie3CSUiomTcIU1ERBfEyoGIiGzH5EBElIG4z4GIiEyCDWkiIkoV2wTHs5WIiCgZKwciIjJxKSsREaUaDj2Hzs5OzJs3D5999hkA4KGHHsKsWbNw22234bbbbsPevXsBAEeOHMHChQsxe/ZsrF+/HpFIxMqwiIgyXsZWDgcPHsSSJUvQ2NhoXjt06BBeffVV1NbWora2FjNnzgQArF27Fhs2bMCePXsghMD27dutCouIaFjI2Mph+/btqKqqgt/vBwD09PTg5MmTqKioQFlZGbZu3Qpd19HU1IRgMIjJkycDABYuXIi6ujqrwiIiymh2ncrqteqJN27cmPD47NmzuOmmm1BVVYX8/HyUl5djx44dGDduHHw+n/l7Pp8Pzc3NF/16BQUjBx0zAGRne2Nx5A/J8w0FGWOSHccqPRyn9Mg0Tkrsr16vamlcliWHZMXFxXjuuefMx8uXL8ebb76JsWPHQlEU87oQIuFxulpaOqEPwR23w+FovyMQ6Bj0cw2VcDiC7GyvVDHJzOfL51ilgeOUHtnGyeg1hMPaoOJSVaXfL9W2rVY6evQo9uzZYz4WQsDr9aKwsBCBQMC8fvbsWXMqioiIPp8Um+BWrVqF+vr6Qb2QEAJPPvkk2tra0Nvbi9dffx0zZ85EUVERcnJy0NDQAACora1FaWnpoF6LiGi46tvnYO3rpJUcZs6cieeffx6zZ8/Gtm3bcP78+Yt+oZKSEqxcuRJLlizB3LlzcfXVV2PevHkAgJqaGlRXV2POnDno7u7GihUrLvr5h0p3pBuBnhbHXp+IKB1WL2VNq+cwf/58zJ8/Hx9//DH+/Oc/Y9GiRZgyZQqWL1+Oa6+9tt8/u2/fPvPvly1bhmXLlqX8TklJCXbs2HGRoVujqfM0AEAXOlSFewSJSDKybYLTdR0nTpxAY2MjNE1DQUEBHnnkEWzdutXK+BwT0sJOh5DA6m8JRJQZ7JpWSqty2LJlC9544w0UFxdj6dKleOaZZ5CVlYXu7m7MmDEDq1evtjZKB4S0EEZ4c50OAwAgQl0InvoUevsZqJexWU9Ekkwrtba24sUXX0RJSUnC9by8PPz617+2JDCnabrmdAgm0XUOAKC1fMrkQEQArN8El9a0kqZpKYnBqBZuueWWoY9KApqQJznA2PchdGfjICJpOFo5VFVVobm5GQ0NDWhtbTWvRyIRfPrpp5YG5rSIRJVDH/YdiNzOSAqO9hwWLVqEY8eO4ejRo5g9e7Z53ePxmGchDVeajN/SdQljIiJHOHq20sSJEzFx4kTcfPPNuOqqqywORS6akPDYcCmrGSJygqPTSj//+c/xzDPP4M477/zcn7/99tuWBCUDGaeVhC5hwiIiWxkpYQiOkutXv8nhrrvuAgA8/PDD1kYhIV2mhrRBY3Igohgnb/ZzzTXXAABuuOEGjB49GjfccAO6u7vxz3/+E1dffbWlgTlNxsoBrByISKYd0hs2bMCLL76Ijz/+GJWVlfjss89QUVFhbWQOk2opa4zQ5IuJiOwT32eQ4lTWQ4cO4ZFHHsHevXuxYMECVFdXo6mpydLAnCbTJjgTKwciV4uvFqSoHIQQUFUV//jHP3DTTTcBAILBoKWBOS0iYeUArdfpCIjIQXp85eBkz8Hwla98BXfddRc+++wz3HDDDXjggQcwfvx4SwNzmi7TPgfjTSBjNUNEjpDiHtLV1dXYu3cvrr/+emRlZWHq1Km4/fbbLQ7NWVIlh9jbQDA5ELlaQs9BhmmlvLw8TJ06Fe3t7Th8+DCuvfZaHD9+3NrIHCZVcmDlQERI3Nugy3Aq6zPPPIM//OEPKCgoMK8pioJ3333XssCcJlVyMNeuyRQTEdktoc8gw/0camtr8c4777jqCA2pkoO5JZKrlYjcLHG1kgQN6dGjR7sqMQDWl2wXhz0HIkpMCI4en2GYNm0aNm/ejO9+97vIze27O9qECRMsC8wJ8dWCVMdnsOdARLD30P60ksMbb7wBAKirqzOvDceeQ0JykOjeCYLJgYhg77RSWslh3759lgYhi/h7OOhS3TvBSA4yxUREdjOnuxXrp5XS6jl0dXXhsccew49+9COcP38eGzZsQFdXl7WROSB+Kkmqm/0Yd35iQ5rI1eJyA6yeZEorOTzxxBPIz89HS0sLcnJy0NnZiQ0bNlgamBM0PX5aSaLkAE4rEVHfVJKiKHJUDkeOHMGaNWvg9XoxYsQI1NTU4MiRI9ZG5oCEaSUJKwcmByJ3k24pq6om/pqmaSnXhoP4aSUmByKSjVk5AHJsgvvWt76Fp556CsFgEH//+9/x6quv4sYbb7Q2MgfIWDlE3wzc50BEcZWDLA3pX/7yl8jLy0N+fj6efvpplJSUYN26ddZG5gBNxsohfr+FLDERkSPiKwfHl7Lu3bsX27Ztw9GjR5Gbm4vx48fjuuuuQ05OjqWBOUGXsHJImEriaiUiVzM/lRTF8tfqNzns3r0bW7ZswerVq1FSUgJFUfDvf/8bGzduRCgUwqxZsywP0E7xd3+TJjlocQmB00pErhZfOTh6KuvLL7+Ml156CV/60pfMa2PHjsWkSZNQUVEx/JKDhNNK8X0G9hyI3E2a24R2dXUlJAbDV7/6VYRCIcuCcoqU00qsHIgopm+fQ+JjK/SbHDwezwV/ZnUzxAmJq5Uk+eczE4LC4zOIXC75Y8nKT6nht1lhEKQ8PsNoQisqG9JELhe/Qzr+sRX67TkcPXoU1113Xcp1IQTC4fCAT97Z2YnFixfjd7/7Hb785S+jvr4e1dXVCIVCuPXWW7FmzRoA0R3Y69evR1dXF6ZOnYpHH30UXm9aWzCGVPzxGUKS4zOM85QUVWHPgcjlUioHC0uHfj+B9+7de8lPfPDgQVRWVqKxsREAEAwGUVFRgVdeeQWjR49GeXk59u/fj+nTp2Pt2rV44oknMHnyZFRUVGD79u1YunTpJb/2pdJkrBy0WEyKyp4Dkcvp8TukYW1y6HdaqaioqN//92f79u2oqqqC3+8HAHz44YcYM2YMiouL4fV6UVZWhrq6OjQ1NSEYDGLy5MkAgIULFybcN8JOMu6QNqaSFCYHIteL3yEdu2LZa1k2d7Nx48aEx2fOnIHP5zMf+/1+NDc3p1z3+Xxobm6+6NcrKBh56cHGjAxmA4iOe1aWCp8vf9DPOVg9PdGYoEaTgwwxZQKOU3o4TumRZZzaQtEviGqs5/DFgpHIzbbmY9y2iX1d180mChDtWyiKcsHrF6ulpRP6IA8bOdcWvUeFoijoCYURCHQM6vmGQqQlFoOiAkLHmTPtlzQ+buLz5Uvx7052HKf0yDROra3Rzyjjky4Q6Ljk5KCqSr9fqm1brVRYWIhAIGA+DgQC8Pv9KdfPnj1rTkXZzdghrUCRcFrJWNjMqSUit0q82Y+DPYehNGnSJHzyySc4ceIENE3Drl27UFpaiqKiIuTk5KChoQEAUFtbi9LSUrvCSmAkBEVRpUkOIn4pK8C+A5GLCRsb0rZNK+Xk5GDTpk1YtWoVQqEQpk+fjjlz5gAAampqUFlZic7OTkyYMAErVqywK6wExmoluSqHWDJQmRyI3K6vIW3cKDQDG9KGffv2mX8/bdo0vPXWWym/U1JSgh07dlgdyoD6KgeJkoMWt1oJ0fOV2HEgcidplrK6jbGUVZWxcuC0EpHriaSmg2NnK7mNOa2kKNKcrSSMg/eMhjSTA5FrDcuGdCbQY8dnRHsOknwIm8dnsHIgcruUysHC12JyiJNQOVh99+50cVqJiGKMyW4F1h+8x+QQJ361kixnKyUvZeXhe0TuZedSViaHOJrQocT+JyRJDkaloLDnQOR6el/pAICVg200EVsmqkh0KqvxbuC0EpHrsXJwiK7rgKJIuAlO6Xs3MDkQuVbKJjhWDvYwKwfIc2S30CPRN4LxZmByIHKtlMrBwtdicogTTQ7RdQCyJIeUSUYmByLXMg6eVthzsJcWm1aCTMdn6Fp0qosNaSLXS04G7DnYJLpaCXJVDiLWc2DlQOR6fS0HJeGxFZgc4uixaSVAnuMzjMqhrwHF5EDkVqmVA6eVbKGJ2BQOIM3xGSkNaFYORK6VfCrrIG9+2S8mhzh9q5Vk6jnoiC5l5bQSkdv1LWVNvjD0mBzimDukFUh0tlIkOr/I5EDkerpuVA7GNLN1r8XkEEePfUs3NsFZOZ+XNtEXU/QhkwORW4nkpawWvhaTQxxNaLFBj468DFNLQtdii5VYORC5HRvSDonfBBd97HxyMI/PSHhMRG6UdHoGp5XsYvQcjJGXYsVS0lJWJgci99KTOtLCwoklJoc4emwKR6rzlcyb/bDnQOR2KT0HVg72iL+fg/HYacKoHLhDmsj1knsMVm7W9Vr2zBko/n4OxmPHGauVzGmliKPhEJFzkm4hbelyJVYOcTTjkDuJVivBOLIbAFQ1liyIyI3MHdIK9znYKmVaSZfggzg+BtXDngORiyUnAyunlZgc4uhxN/sxHjtNxC9lVT3sORC5mLlDWhngF4cAk0McTegJ906QoSFtLmUFoKheJgciF0teuspNcDbRzCO7jceSJAeD6mFDmsjF+payGn1R616LySFO381+JNoEF6tmAACKCiFDH4SIHJF8D2meymoTPTa/L920EnsORIS4SsFY2W7hazE5xNGE3refAHIsZRVCN98ICpMDkaslVw7sOdhAFzoEhHzTSsbNfgCADWkiV+vLBdznYBtjCinhVFYZ5vdFYkNasCFN5Fpm5WCc4mBhR5rJIUYzD7iDOfKyHJ9hTnVxWonI1XSBxMNBLUwOjpyttHz5crS2tsLrjb78Y489hq6uLlRXVyMUCuHWW2/FmjVrbI0pIqLfyOMrBxl6Dok7pHl8BpGbCSESz38bTslBCIHGxkb87W9/M5NDMBjEnDlz8Morr2D06NEoLy/H/v37MX36dNviisSmaxRFnlNZhRDmwXsAG9JEbpd8tpKmWfcZZXtyOH78OADgJz/5Cc6fP48f/OAH+MY3voExY8aguLgYAFBWVoa6ujpbk4MxrSTVzX6SD29XvRCRkHPxEJGjkk9lHVaVQ3t7O6ZNm4aHH34Yvb29WLFiBe688074fD7zd/x+P5qbmy/qeQsKRg4qrt72LgBAlteDLK8HAJA3Mhs+X/6gnncwRKQXnQBUT7Q1lJ2bDb076GhMmYJjlB6OU3pkGafc3CwoioKsLA8EgLwv5FgWm+3JYcqUKZgyZYr5eNGiRdi6dSuuv/5685oQImG/QTpaWjoH1Zw509kOANA0gUhsa8n59m4EAh2X/JyDZVQJuhb95+qNAHo47GhMmcDny+cYpYHjlB6Zxqm7OwwAiER0eACcb+u55NhUVen3S7Xtq5X+9a9/4YMPPjAfCyFQVFSEQCBgXgsEAvD7/bbGZfYcEhrSDk8rGc1oIyBFTWxQE5Gr6MlLWS3sOdieHDo6OrB582aEQiF0dnZi586d+MUvfoFPPvkEJ06cgKZp2LVrF0pLS22NK2L0HOJPZXX6g9hsiPctZeX9HIjcS4jYavvY42HVc5gxYwYOHjyI22+/HbquY+nSpZgyZQo2bdqEVatWIRQKYfr06ZgzZ46tccVXDjB3SDu8Wil+7wXAfQ5ELidi2aHvC+wwSg4AcP/99+P+++9PuDZt2jS89dZbToQDAIiIuMohds3xaaWUyoHHZxC5mUg6eI87pG3QVzlAnlNZ9cTkoPB+DkSulnzw3rDqOcgqcVopyvHkIFKnldhzIHIv3dghDUBVFFYOdtASppUkOVspuYZUuVqJyM2iZytFPw88HiYHWyRXDgoUxxvSKf0FNqSJXE3EVw6qYunBe0wOMfFnKwGAR1EdTw7mLUEVo+fAhjSRm8Xfv8GrKtA0JgfLmauVYnlZVdS+Y7ydkpycVA8gNEvv/kRE8or/b9+jKpbuxWJyiEmuHFTF43jl0Hejn7j7OSRcJyI3EbH7OQCAx6Miwmkl6yX3HDyq6vxqpeTjM1Q18ToRuUryaiX2HGyg6VrCMlZVUSVYrfQ5+xwA9h2IXCp+RpmrlWzSKyLwqn0bxj2Kx/HkID5vhzTA5EDkUiJuXsmjKtwEZ4ew1oscT7b52Kt6nW9If97ZSgAEd0kTuZIu+k5k9agqKwc7hLUwstQs87FX9aLX6Q9hPfVU1uh1Vg5EbhTfY+C0kk3CWjihcshSvYjovQ5GBLPnYHxTYM+ByN3iG9LRpaxMDpYL673I9vRVDtHk4HRDOmkpq8LVSkRupuuibykrew72iE4rxfUcFC96na4ckr8VxBrSPHyPyJ1YOTggpSHtcb7nIMxTWZN7DmxIE7lRtOdgHLzHhrQtQno4aVopy9wY5xjucyCiOLou4lYrsXKwRa8WRnZ85aB4nE8Oyb0FJgciV4vPBR4evGePoBZKSA5Zapbj00p9+xwSp5XYcyByJy2lcmBD2lIRPYKu3m5cljXSvBbtOTjdkP6c+zkAqae1EpEr6II9B1t1hDsBAJflXGZek2Epq7kTWknuObAhTeRGQufBe7ZqC7cDAC7PzjeveRUJNsHpSfscYg1zEXE4LiJyhCbi9jlwh7T12kIdAIDLkysHoTl7TwctqXLIyo0+7g06FBAROUmPqxy83ARnvfZY5XBZXOVgnLPk5OF7IrkhHUsOgsmByJUSd0iz52C5tlAHFCgJycEbm993dMWS+dqxyiGbyYHIzaI7pI2GtMI7wVmtPdyBL2TlwWM0fAF4Y5WDo8lBiwBKX0zwZEerCCYHIlfS9fgbQ7Ihbbn2cHtCvwGAuVs6rIWdCAlAbLWSpy85KIoCeHNZORC5lC6SD95jcrBUW6gjYUoJAHK90SmcoBZyIqQoXeu7+1uMkj0CItztUEBE5KSE+zmoCnQhYneHG3pMDohOK12enVg55HpyAADBiIPf0rUIFE9Schj5RYjOVocCIiInRXdI922CM65ZwfXJQRc62sMduCwnsXIYYVYODiYHPdK3KzpGzb8SekfAoYCIyEkRTU84shtgcrBMV283dKH3Uzk4N60kImHAm5NwTc33QXS2QGjcCEfkNhFNTzhbCYBlfQfXJ4fzoTYAwBVJDelcCSoH0RuEkpWUHEYVAUJAbzvtUFRE5AQhBCJa3LSSWTlYsxHO9cmhNXgeADAq94qE6zkSVA7oDfXtio5RRxUBAPTWJiciIiKHRLTEe8qz52CxcxdMDtlQoKDHwYa0iIRSp5WuKAQUFfo5JgciN+mNRJOAmlQ5WLXXQark8Pbbb+P73/8+Zs2ahddee82W1zwXOg+v6sXIrC8kXFcUBVfkXG5OOzlBhHtSKgfFkwX1Mj/0cycdioqInNCbXDnEkoNVu6S9A/+KPZqbm7Flyxa88cYbyM7OxuLFi3HjjTfi61//uqWve6qrGf4RV0JVUvPkF3NHoSXo3LJR0X0eSvHElOvqqCJorByIXKWrJ7oIRY0lhZys6ErGYMiaUxykSQ719fW46aabcMUVVwAAZs+ejbq6Otx3331p/XljwC5Gd28P2sJtmOi72vzzfr/ffL5xX/wqDp09Ag0aslR7h0rvbIU3Lx9Z/v+D3+9HVpbHjDGr6BvAuf9B0VJ7EhR1Ke8HN+I4pUeGceoJR+AfNQKXjS7EiBwvRl/5BfhHjUBnMHJJ8Q30ZxRh1fa6i/T73/8e3d3dWLNmDQDgT3/6Ez788EM8/vjjDkdGROQ+0vQcdF03l2gB0WVb8Y+JiMg+0iSHwsJCBAJ9O38DgYA5xUNERPaSJjl8+9vfxgcffIDW1lb09PTgnXfeQWlpqdNhERG5kjQN6auuugpr1qzBihUr0Nvbi0WLFuHaa691OiwiIleSpiFNRETykGZaiYiI5MHkQEREKZgciIgoBZMDERGlYHIgIqIUTA6SGehk2r/+9a+47bbbMH/+fNxzzz1oa3Pu1FgnpXuC73vvvYfvfOc7NkYmn4HG6vjx41i+fDnmz5+Pn/70p3xPXWCcDh8+jDvuuAPz589HeXk52tvbHYjSRoKkcfr0aTFjxgxx7tw50dXVJcrKysSxY8fMn3d0dIibb75ZnD59WgghxNNPPy0ef/xxp8J1zEDjZAgEAmLOnDlixowZDkQph4HGStd1MWvWLLF//34hhBBPPfWU2Lx5s1PhOiad99SSJUvEe++9J4QQorq6WvzmN79xIlTbsHKQSPzJtHl5eebJtIbe3l5UVVXhqquuAgCMHz8ep06dcipcxww0TobKysq0T/UdrgYaq8OHDyMvL888jeDuu+/GsmXLnArXMem8p3RdR1dXFwCgp6cHubnD+0RkJgeJnDlzBj6fz3zs9/vR3NxsPh41ahRmzpwJAAgGg3jhhRfwve99z/Y4nTbQOAHAyy+/jG9+85uYNGmS3eFJZaCx+t///ocrr7wSFRUVWLBgAaqqqpCXl+dEqI5K5z314IMPorKyErfccgvq6+uxePFiu8O0FZODRNI9mbajowMrV65ESUkJFixYYGeIUhhonD766CO88847uOeee5wITyoDjVUkEsGBAwewZMkS7Ny5E8XFxdi0aZMToTpqoHEKBoNYv349XnrpJbz//vtYunQpfvWrXzkRqm2YHCSSzsm0Z86cwdKlSzF+/Hhs3LjR7hClMNA41dXVIRAI4I477sDKlSvNMXOjgcbK5/NhzJgxmDgxesfBefPm4cMPP7Q9TqcNNE4fffQRcnJyzPPefvjDH+LAgQO2x2knJgeJDHQyraZpuPvuu3Hrrbdi/fr1rr3fxUDjtHr1auzZswe1tbV44YUX4Pf78cc//tHBiJ0z0FhNmTIFra2t+O9//wsA2LdvHyZMmOBUuI4ZaJzGjBmD06dP4/jx4wCAd99910yow5U0p7LShU+mveuuu7B69WqcPn0a//nPf6BpGvbs2QMAuOaaa1xXQQw0TsP9P9qLkc5YPffcc6isrERPTw8KCwuxefNmp8O2XTrjVF1djfvvvx9CCBQUFODJJ590OmxL8VRWIiJKwWklIiJKweRAREQpmByIiCgFkwMREaVgciAiohRMDkRElILJgYiIUvx/iZJ2xkV/2VgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "basis_function = Polynomial(degree=2)\n",
    "parameters = np.zeros([3, 50])\n",
    "\n",
    "for i in range(50):\n",
    "    x_train, x_valid, y_train, y_valid = get_siso_data(n=3000,\n",
    "                                                       colored_noise=True,\n",
    "                                                       train_percentage=90)\n",
    "\n",
    "    model = FROLS(\n",
    "        order_selection=False,\n",
    "        n_terms=3,\n",
    "        extended_least_squares=True,\n",
    "        ylag=2, xlag=2, elag=2,\n",
    "        info_criteria='aic',\n",
    "        estimator='least_squares',\n",
    "        basis_function=basis_function\n",
    "    )\n",
    "    \n",
    "    model.fit(X=x_train, y=y_train)\n",
    "    parameters[:, i] = list(model.theta)\n",
    "\n",
    "sns.set()\n",
    "pal = sns.cubehelix_palette(3, rot=-.5, dark=.3)\n",
    "\n",
    "ax = sns.kdeplot(parameters.T[:, 0])\n",
    "ax = sns.kdeplot(parameters.T[:, 1])\n",
    "ax = sns.kdeplot(parameters.T[:, 2])\n",
    "# plotting a vertical line where the real values must lie\n",
    "ax = plt.axvline(x=0.1, c='k')\n",
    "ax = plt.axvline(x=0.2, c='k')\n",
    "ax = plt.axvline(x=0.9, c='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Now we have an unbiased estimation of the parameters! \n",
    "\n",
    "## Note\n",
    "\n",
    "Note: The Extended Least Squares is an iterative algorithm. In SysIdentpy we fixed 30 iterations because it is known from literature that the algorithm converges quickly (about 10 or 20 iterations)."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "0e65fe37feb8ff9f7778552a28949e943d61f86c936833305e2c18cda5b438ac"
  },
  "kernelspec": {
   "display_name": "Python 3.8.11 64-bit ('rd': conda)",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}