regression_multiple_to_multiple.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression \n",
    "multiple predictors to multiple predictants\n",
    "\n",
    "$$\\begin{bmatrix} y_1 & y_2 & y_3 \\end{bmatrix} \n",
    "= \\begin{bmatrix} x_1 & x_2 & 1 \\end{bmatrix} \n",
    "\\begin{bmatrix} \n",
    "    \\beta_{1, 1} & \\beta_{2, 1} & \\beta_{3, 1} \\\\ \n",
    "    \\beta_{1, 2} & \\beta_{2, 2} & \\beta_{3, 2} \\\\ \n",
    "    \\beta_{1, 0} & \\beta_{2, 0} & \\beta_{3, 0} \\\\ \n",
    "\\end{bmatrix} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import signal, linalg\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x.shape = (1008, 2)\n",
      "y.shape = (1008, 3)\n"
     ]
    }
   ],
   "source": [
    "# Read Nino indices\n",
    "f = 'data/era5_nino.csv'\n",
    "nino = pd.read_csv(f, index_col=0, parse_dates=True)\n",
    "\n",
    "# Read T2m data\n",
    "f = 'data/era5_monthly_t2m_points.csv'\n",
    "df = pd.read_csv(f, index_col=0, parse_dates=True)\n",
    "\n",
    "# Estimate anomaly\n",
    "dfa = df.groupby(df.index.month).transform(lambda x: x-x.mean())\n",
    "\n",
    "# Detrend\n",
    "data = signal.detrend(dfa.values, axis=-1)\n",
    "\n",
    "# Define predictor and predictand\n",
    "x = nino[['nino3', 'nino4']].values\n",
    "y = data[:, :3]\n",
    "n = x.shape[0]\n",
    "\n",
    "print(f'{x.shape = :}')\n",
    "print(f'{y.shape = :}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$y\n",
    "= \\begin{bmatrix} y_1 & y_2 & y_3 \\end{bmatrix} \n",
    "= Ap\n",
    "= \\begin{bmatrix} x_1 & x_2 & 1 \\end{bmatrix} \n",
    "\\begin{bmatrix} \n",
    "    \\beta_{1, 1} & \\beta_{2, 1} & \\beta_{3, 1} \\\\ \n",
    "    \\beta_{1, 2} & \\beta_{2, 2} & \\beta_{3, 2} \\\\ \n",
    "    \\beta_{1, 0} & \\beta_{2, 0} & \\beta_{3, 0} \\\\ \n",
    "\\end{bmatrix} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A.shape = (1008, 3)\n",
      "[[0.77482826 0.49944904 1.        ]\n",
      " [0.91584975 0.6616926  1.        ]\n",
      " [0.96427315 0.78748566 1.        ]\n",
      " [0.48696518 0.6775541  1.        ]\n",
      " [1.1274626  0.87392026 1.        ]]\n"
     ]
    }
   ],
   "source": [
    "A = np.hstack([x, np.ones((n, 1))])\n",
    "print(f'{A.shape = :}')\n",
    "print(A[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p.shape = (3, 3)\n",
      "p = [[-6.51854590e-02  1.68140215e-01 -1.56988096e-01]\n",
      " [ 1.28367018e-01 -2.26603304e-01  4.44847279e-02]\n",
      " [ 8.81011244e-08 -1.55595957e-07  3.07159735e-08]]\n",
      "Sums of squared residuals = [ 695.95524571 1568.9124786  1554.14069974]\n",
      "Rank of matrix A (number of linearly independent columns): 3\n",
      "Singular values of A: [31.74901573 29.80843068  9.99211823]\n"
     ]
    }
   ],
   "source": [
    "p, res, rnk, s = linalg.lstsq(A, y)\n",
    "\n",
    "print(f'{p.shape = :}')\n",
    "print(f'{p = :}')\n",
    "print(f'Sums of squared residuals = {res}')\n",
    "print(f'Rank of matrix A (number of linearly independent columns): {rnk}')\n",
    "print(f'Singular values of A: {s}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y_res.shape = (1008, 3)\n",
      "Sums of squared residuals = [ 695.95524571 1568.9124786  1554.14069974]\n"
     ]
    }
   ],
   "source": [
    "# Reconstruction\n",
    "y_recon = A @ p\n",
    "\n",
    "# Residual\n",
    "y_res = y - y_recon\n",
    "print(f'{y_res.shape = :}')\n",
    "\n",
    "# sum of squared residuals\n",
    "print(f'Sums of squared residuals = {(y_res**2).sum(axis=0)}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "basicf9",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}