{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "62fa0297",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "\n",
    "@author: Gabriel Turinici 2025\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "This program uses historical data in the format in :\n",
    "    https://turinici.com/wp-content/uploads/cours/common/close_cac40_historical.csv\n",
    "\n",
    "It can also be downloaded form yahoo finance in daily\n",
    "prices (at least the \"close\") if possible at lest 5 years\n",
    "\n",
    "Idea: use yahoo e.g., yfinance package\n",
    "\"pip install yfinance\"\n",
    "\n",
    "\n",
    "Then the code does :\n",
    "\n",
    "1'/ order by increasing date\n",
    "2/ plot price histogram and returns  (with \"log\" and/or \"actuarial\")\n",
    "3/ test normality of : prices, log returns, actuarial returns\n",
    "    for instance can use scipy.stats.normaltest\n",
    "4/ shows the random versus optimal results\n",
    "5/ show what a backtest would obtain : compute market portfolio from T-6Y,T-1Y data\n",
    "(last 5 years of data before this last year)\n",
    "6/ compute performance of estimated market portfolio on the T-1Y,T data (last year data)\n",
    "7/ compare with random choices\n",
    "\n",
    "\"\"\"\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import kstest, normaltest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "3202bb04",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 166864 characters]",
      "text/plain": [
       "<Figure size 1500x1500 with 42 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# from scipy.special import softmax\n",
    "\n",
    "# take csv from www course :we suppose it is available locally\n",
    "data = pd.read_csv(\n",
    "    \"https://turinici.com/wp-content/uploads/cours/common/close_cac40_historical.csv\",\n",
    "    sep=\";\",\n",
    "    index_col=\"Date\",\n",
    ")\n",
    "data.head()\n",
    "# order by increasing date, keep variable 'data'\n",
    "data.sort_index(ascending=True, inplace=True)\n",
    "# data.sort(by='Date',inplace=True)\n",
    "data.head()\n",
    "data.tail()\n",
    "\n",
    "# plot histogram of prices\n",
    "_ = data.hist(bins=30, figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "f1d15e71",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test pval= 2.0507269291153243e-147 res= not normal\n",
      "Test pval= 1.2341507304168457e-145 res= not normal\n",
      "Test pval= 1.5109663099066784e-82 res= not normal\n",
      "Test pval= 5.193407329538019e-261 res= not normal\n",
      "Test pval= 9.263758145504221e-171 res= not normal\n",
      "Test pval= 1.5679377988996054e-95 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 3.431480107389548e-229 res= not normal\n",
      "Test pval= 7.209347413730216e-130 res= not normal\n",
      "Test pval= 5.208576684311371e-248 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 6.226583000372112e-32 res= not normal\n",
      "Test pval= 2.8139918420327365e-122 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 3.2023586082752444e-62 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 2.3734732069569748e-122 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 3.658142667687734e-131 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 9.297487741263689e-210 res= not normal\n",
      "Test pval= 4.5583316571649815e-156 res= not normal\n",
      "Test pval= 7.965502395433165e-197 res= not normal\n",
      "Test pval= 4.631868468721008e-226 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "Test pval= 5.162726389052399e-169 res= not normal\n",
      "Test pval= 2.1068459098663725e-95 res= not normal\n",
      "Test pval= 1.724502624699572e-235 res= not normal\n",
      "Test pval= 6.864779304946272e-62 res= not normal\n",
      "Test pval= 1.701124981576563e-154 res= not normal\n",
      "Test pval= 4.3391395462980045e-304 res= not normal\n",
      "Test pval= 0.0 res= not normal\n",
      "no. of normal =  0\n",
      "no. of not normal =  40\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0, 40)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def normality_test(data, kolmogorov_smirnov=False, level=0.01, print_results=True):\n",
    "    \"\"\"\n",
    "    Tests normality of each column of dataframe \"data\".\n",
    "    Inputs:\n",
    "        kolmogorov_smirnov= false: use \"normaltest\", otherswise use kstest, both from scipy.stats\n",
    "        level = p-value threshold level for the conclusions\n",
    "\n",
    "    Outputs:  the number of yes/no in the results\n",
    "    \"\"\"\n",
    "    pvalues = []\n",
    "    for cols in data.keys():\n",
    "        pv = (\n",
    "            kstest(data[cols].dropna(), \"norm\").pvalue\n",
    "            if kolmogorov_smirnov\n",
    "            else normaltest(data[cols].dropna()).pvalue\n",
    "        )\n",
    "        pvalues.append(pv)\n",
    "        if pv < level:\n",
    "            res = \"not normal\"\n",
    "        else:\n",
    "            res = \"normal\"\n",
    "        print(\"Test pval=\", pv, \"res=\", res)\n",
    "    normalok = sum([1 for pv in pvalues if pv >= level])\n",
    "    normalnotok = sum([1 for pv in pvalues if pv < level])\n",
    "    if print_results:\n",
    "        print(\"no. of normal = \", normalok)\n",
    "        print(\"no. of not normal = \", normalnotok)\n",
    "    return normalok, normalnotok\n",
    "\n",
    "\n",
    "normality_test(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "fba88600",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/arthurdanjou/Workspace/studies/.venv/lib/python3.13/site-packages/pandas/core/internals/blocks.py:393: RuntimeWarning: invalid value encountered in log\n",
      "  result = func(self.values, **kwargs)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test pval= 0.025428757930302828 res= normal\n",
      "Test pval= 0.01412845580561908 res= normal\n",
      "Test pval= 0.2844927942646553 res= normal\n",
      "Test pval= 0.2807546954064213 res= normal\n",
      "Test pval= 0.009428521835360353 res= not normal\n",
      "Test pval= 0.14982216104037197 res= normal\n",
      "Test pval= 0.07490385366598394 res= normal\n",
      "Test pval= 0.1693491210857498 res= normal\n",
      "Test pval= 0.00622835412057237 res= not normal\n",
      "Test pval= 0.0038601712111598595 res= not normal\n",
      "Test pval= 0.00014234445964435524 res= not normal\n",
      "Test pval= 0.1559244997516793 res= normal\n",
      "Test pval= 5.4492134141853295e-12 res= not normal\n",
      "Test pval= 0.027428045061480538 res= normal\n",
      "Test pval= 0.00011538203750856845 res= not normal\n",
      "Test pval= 0.07350032657006557 res= normal\n",
      "Test pval= 0.10474460687572516 res= normal\n",
      "Test pval= 0.0015506993297124128 res= not normal\n",
      "Test pval= 0.36679129167153646 res= normal\n",
      "Test pval= 0.011263558096640853 res= normal\n",
      "Test pval= 1.4403372640862255e-05 res= not normal\n",
      "Test pval= 0.017985610179734428 res= normal\n",
      "Test pval= 0.008018485213850552 res= not normal\n",
      "Test pval= 0.06774855935762493 res= normal\n",
      "Test pval= 0.057097963780292554 res= normal\n",
      "Test pval= 0.019022081578872616 res= normal\n",
      "Test pval= 0.021072819215708735 res= normal\n",
      "Test pval= 0.5847583599406958 res= normal\n",
      "Test pval= 0.34703248983489265 res= normal\n",
      "Test pval= 0.005743759434418153 res= not normal\n",
      "Test pval= 0.05305732031736986 res= normal\n",
      "Test pval= 0.0009839608073665932 res= not normal\n",
      "Test pval= 0.07283923280212494 res= normal\n",
      "Test pval= 2.232637543452161e-06 res= not normal\n",
      "Test pval= 0.012943908206025189 res= normal\n",
      "Test pval= 0.05971194449910225 res= normal\n",
      "Test pval= 0.002870485543548689 res= not normal\n",
      "Test pval= 0.042672987739536226 res= normal\n",
      "Test pval= 0.029170249033746622 res= normal\n",
      "Test pval= 0.002812593222739167 res= not normal\n",
      "no. of normal =  26\n",
      "no. of not normal =  14\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(26, 14)"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 156424 characters]",
      "text/plain": [
       "<Figure size 1500x1500 with 42 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use 'data' to compute returns\n",
    "# returns= data.pct_change() #actuarial\n",
    "returns = np.log(data / data.shift(1))\n",
    "\n",
    "nr_days = 252\n",
    "\n",
    "# print('annualize returns')\n",
    "# returns *= nr_days\n",
    "\n",
    "_ = returns.hist(bins=int(np.sqrt(returns.shape[0])), figsize=(15, 15))\n",
    "\n",
    "normality_test(returns.tail(25 * 3))   # test last 3 months"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "5b14395c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "normality tests for increments, not returns!!\n",
      "Test pval= 0.08683895238128052 res= normal\n",
      "Test pval= 0.021035773553864635 res= normal\n",
      "Test pval= 0.38286247315996325 res= normal\n",
      "Test pval= 0.3710853667901102 res= normal\n",
      "Test pval= 0.0395519075986371 res= normal\n",
      "Test pval= 0.26910192261508514 res= normal\n",
      "Test pval= 0.11094121314729449 res= normal\n",
      "Test pval= 0.3180825721276856 res= normal\n",
      "Test pval= 0.011858172793204136 res= normal\n",
      "Test pval= 0.008952882580561775 res= not normal\n",
      "Test pval= 0.0008467555154087856 res= not normal\n",
      "Test pval= 0.2098816393496899 res= normal\n",
      "Test pval= 5.524153122423706e-14 res= not normal\n",
      "Test pval= 0.15309417383852073 res= normal\n",
      "Test pval= 0.0010024995182926846 res= not normal\n",
      "Test pval= 0.3549281579113689 res= normal\n",
      "Test pval= 0.3129404747848908 res= normal\n",
      "Test pval= 0.008997473704571489 res= not normal\n",
      "Test pval= 0.6147399643675131 res= normal\n",
      "Test pval= 0.09700094185831638 res= normal\n",
      "Test pval= 0.00118024612616608 res= not normal\n",
      "Test pval= 0.010093733928468086 res= normal\n",
      "Test pval= 0.03923585485509744 res= normal\n",
      "Test pval= 0.15164264770844182 res= normal\n",
      "Test pval= 0.2480354891922463 res= normal\n",
      "Test pval= 0.0183314268046639 res= normal\n",
      "Test pval= 0.011137481909481228 res= normal\n",
      "Test pval= 0.4300901492943915 res= normal\n",
      "Test pval= 0.5892009761977898 res= normal\n",
      "Test pval= 0.0008843165602394642 res= not normal\n",
      "Test pval= 0.20729225915186456 res= normal\n",
      "Test pval= 0.016863594683601053 res= normal\n",
      "Test pval= 0.14688927380791153 res= normal\n",
      "Test pval= 2.380301465340281e-05 res= not normal\n",
      "Test pval= 0.029281292837861225 res= normal\n",
      "Test pval= 0.20712701484300206 res= normal\n",
      "Test pval= 0.015888165155116662 res= normal\n",
      "Test pval= 0.096035363594102 res= normal\n",
      "Test pval= 0.07771383799032755 res= normal\n",
      "Test pval= 0.004196457641565341 res= not normal\n",
      "no. of normal =  31\n",
      "no. of not normal =  9\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(31, 9)"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 169888 characters]",
      "text/plain": [
       "<Figure size 1500x1500 with 42 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"normality tests for increments, not returns!!\")\n",
    "increments = data - data.shift(1)\n",
    "_ = increments.hist(bins=int(np.sqrt(increments.shape[0])), figsize=(15, 15))\n",
    "\n",
    "normality_test(increments.tail(25 * 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "27153b45",
   "metadata": {},
   "outputs": [],
   "source": [
    "nb = 5  # will work with nb stocks\n",
    "all_returns = returns.copy()  # backup\n",
    "nb_all = all_returns.shape[1]\n",
    "\n",
    "if nb > nb_all:\n",
    "    print(\"too many number of stocks, revert to max\")\n",
    "    nb = nb_all\n",
    "\n",
    "# choose the stock names\n",
    "nb_stocks_names = np.random.choice(all_returns.keys(), nb, replace=False)\n",
    "returns_small = all_returns.loc[:, nb_stocks_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "4e01bf0f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 25900 characters]",
      "text/plain": [
       "<Figure size 600x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mean_returns = returns_small.mean()\n",
    "cov_matrix = returns_small.cov()\n",
    "rdt_list = []\n",
    "std_list = []\n",
    "\n",
    "for _ in range(500):\n",
    "    # sample at random some \"allocation\"\n",
    "    allocation = np.random.dirichlet(np.ones(nb))\n",
    "    rdt_port = allocation @ mean_returns\n",
    "    std_port = np.sqrt(allocation @ cov_matrix @ allocation)\n",
    "    rdt_list.append(rdt_port)\n",
    "    std_list.append(std_port)\n",
    "\n",
    "std_list = np.array(std_list)\n",
    "rdt_list = np.array(rdt_list)\n",
    "\n",
    "inverse_cov = np.linalg.inv(cov_matrix)\n",
    "\n",
    "# compute and draw the efficient frontier on the same graph\n",
    "onesM = np.ones_like(mean_returns)\n",
    "\n",
    "# compute 'a' and 'b' using formulas from the course\n",
    "a = onesM @ inverse_cov @ onesM  # np.sum(inverse_cov)\n",
    "b = onesM @ inverse_cov @ mean_returns\n",
    "\n",
    "# plot the frontier and its symmetric w/r to origin\n",
    "sigmarange = np.linspace(1.0 / np.sqrt(a) + 1.0e-10, 1.1 * np.max(std_list), 47)\n",
    "\n",
    "# compute the return of the optimal portfolio for sigma in sigmarange\n",
    "# will use the \"factor\" auxiliary variable\n",
    "factor = np.sqrt((mean_returns - b / a) @ inverse_cov @ (mean_returns - b / a))\n",
    "optimal_return = b / a + np.sqrt(sigmarange**2 - 1.0 / a) * factor\n",
    "\n",
    "fig = plt.figure(\"perf\", figsize=(6, 3))\n",
    "plt.title(\"efficient frontier\")\n",
    "plt.scatter(std_list * np.sqrt(nr_days), rdt_list * nr_days, s=10)\n",
    "plt.plot(sigmarange * np.sqrt(nr_days), optimal_return * nr_days, \"r-\")\n",
    "plt.xlabel(\"std\")\n",
    "plt.ylabel(\"rdt\")\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "eec4e675",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 11400 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 15772 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Annualized rdt_backtest: -37.97%, std_backtest:330352.97%\n"
     ]
    }
   ],
   "source": [
    "# plot market portfolio\n",
    "r_rate = 0.0  # attention to convert to daily (not annualize!) if not zero !!!\n",
    "market_portf_allocation = inverse_cov @ (mean_returns.values - r_rate * onesM)\n",
    "market_portf_allocation /= np.sum(market_portf_allocation)\n",
    "\n",
    "plt.title(\"market portf\")\n",
    "plt.bar(list(range(1, 1 + nb)), market_portf_allocation)\n",
    "plt.show()\n",
    "plt.close()\n",
    "np.sum(np.abs(market_portf_allocation)) * 100\n",
    "\n",
    "# test on data obtained from T-6Y,T-1Y the rdv/vol on T-1Y,T\n",
    "# and compare with market portfolio on T-1,T. T=last known date\n",
    "\n",
    "# calcul of mean_returns and inverse_cov between T-6Y and T-1Y, 1Y=255 days\n",
    "ti = -6 * 255\n",
    "tf = -255\n",
    "hist_mean_returns = returns_small.iloc[ti:tf, :].mean()\n",
    "hist_cov_matrix = returns_small.iloc[ti:tf, :].cov()\n",
    "hist_inverse_cov = np.linalg.inv(hist_cov_matrix)\n",
    "hist_market_portf_allocation = hist_inverse_cov @ (hist_mean_returns - r_rate * onesM)\n",
    "hist_market_portf_allocation /= np.sum(hist_market_portf_allocation)\n",
    "\n",
    "plt.title(\"estimated market portf\")\n",
    "plt.bar(list(range(1, 1 + nb)), hist_market_portf_allocation)\n",
    "plt.show()\n",
    "plt.close()\n",
    "\n",
    "# compute performance during T-1Y,T\n",
    "c_mean_returns = returns_small.iloc[tf:, :].mean()\n",
    "c_cov_matrix = returns_small.iloc[tf:, :].cov()\n",
    "c_inverse_cov = np.linalg.inv(c_cov_matrix)\n",
    "\n",
    "rdt_backtest = hist_market_portf_allocation @ c_mean_returns\n",
    "std_backtest = (\n",
    "    hist_market_portf_allocation @ c_inverse_cov @ hist_market_portf_allocation\n",
    ")\n",
    "\n",
    "np.sqrt(hist_market_portf_allocation @ c_inverse_cov @ hist_market_portf_allocation)\n",
    "\n",
    "print(\n",
    "    f\"Annualized rdt_backtest: {rdt_backtest * (nr_days) * 100:.2f}%, std_backtest:{std_backtest * np.sqrt(nr_days):.2f}%\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "1ea8bde1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 27152 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nr_tests = 500\n",
    "c_rdt_list = np.zeros((nr_tests,))\n",
    "c_std_list = np.zeros((nr_tests,))\n",
    "\n",
    "for ii in range(nr_tests):\n",
    "    # sample at random some \"allocation\"\n",
    "    allocation = np.random.dirichlet(np.ones(nb))\n",
    "    c_rdt_list[ii] = allocation @ c_mean_returns\n",
    "    c_std_list[ii] = np.sqrt(allocation @ c_cov_matrix @ allocation)\n",
    "\n",
    "# reduce to target volatility (the mean of random ones)\n",
    "target_vol = np.mean(std_list * np.sqrt(nr_days))\n",
    "\n",
    "reduction_factor = target_vol / np.sqrt(nr_days) / std_backtest\n",
    "reduced_hist_allocation = reduction_factor * hist_market_portf_allocation\n",
    "reduced_rdt_backtest = reduced_hist_allocation @ c_mean_returns\n",
    "reduced_std_backtest = np.sqrt(\n",
    "    reduced_hist_allocation @ c_inverse_cov @ reduced_hist_allocation\n",
    ")\n",
    "\n",
    "fig = plt.figure(\"perf last backtest\")\n",
    "plt.scatter(c_std_list * np.sqrt(nr_days), c_rdt_list * nr_days, s=10, label=\"random\")\n",
    "plt.scatter(\n",
    "    reduced_std_backtest * np.sqrt(nr_days),\n",
    "    reduced_rdt_backtest * nr_days,\n",
    "    marker=\"*\",\n",
    "    label=\"backtest\",\n",
    ")\n",
    "plt.xlabel(\"std\")\n",
    "plt.ylabel(\"rdt\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.close()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "studies",
   "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.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
