{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b1cbed9b",
   "metadata": {},
   "source": [
    "<h1> Jupyter-Notebook zu Copulas </h1>\n",
    "(Lisa Kusch, Nov 2022)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1343856f",
   "metadata": {},
   "source": [
    "## 1 Grundlagen: Wahrscheinlichkeitsdichte und Wahrscheinlichkeitsverteilung\n",
    "\n",
    "Die Dichtefunktion der Normalverteilung ist gegeben durch\n",
    "\n",
    "$f(x\\mid \\mu ,\\sigma)={\\frac {1}{\\sqrt {2\\pi \\sigma ^{2}}}}\\operatorname {exp} \\left(-{\\frac {(x-\\mu )^{2}}{2\\sigma ^{2}}}\\right)\\quad$ für $\\quad -\\infty <x<\\infty $.\n",
    "\n",
    "Die zugehörige Verteilungsfunktion ist\n",
    "\n",
    "$F(x\\mid \\mu, \\sigma)= {\\frac {1}{2}}\\left(1+\\operatorname {erf} \\left({\\frac {x-\\mu}{\\sqrt {2}\\sigma}}\\right)\\right)$,\n",
    "\n",
    "wobei $\\operatorname{erf}$ die Fehlerfunktion ist.\n",
    "\n",
    "Im Folgenden sollen diese Funktionen geplottet werden. Dazu wird das Modul `matplotlib` verwendet. Die oben beschriebenen Funktionen stehen bereits in `scipy.stats` zur Verfügung. Diese werden zunächst importiert."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abc7ed3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Benötigte Imports\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "\n",
    "# Seed wird festgesetzt (Reproduzierbarkeit)\n",
    "np.random.seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2827351f",
   "metadata": {},
   "source": [
    "Nun wird die Funktion im Bereich $[-4,4]$ geplottet. Dabei ist $\\mu = 0$ und $\\sigma = 1$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc99cef4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotte Normalverteilung mit Erwartungswert mu und Standardabweichung sigma\n",
    "mu=0.0\n",
    "sigma=1.0\n",
    "\n",
    "#Plotte zwischen -4 und 4 mit Auflösung von 0.005 (Schrittweite)\n",
    "x = np.arange(-4, 4, 0.005)\n",
    "\n",
    "plt.plot(x, norm.pdf(x,mu,sigma), c='black')\n",
    "plt.xlabel('x', fontsize=13)\n",
    "plt.ylabel('f(x)', fontsize=13)\n",
    "plt.show()\n",
    "\n",
    "plt.plot(x, norm.cdf(x,mu,sigma), c='black')\n",
    "plt.xlabel('x', fontsize=13)\n",
    "plt.ylabel('F(x)', fontsize=13)\n",
    "plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41004b33",
   "metadata": {},
   "source": [
    "### Aufgabe 1.1:\n",
    "Probieren Sie verschiedene Werte von $\\mu$ und $\\sigma$ aus. Wie verändert sich die Funktion? Was bedeuten die beiden Parameter?\n",
    "\n",
    "### Aufgabe 1.2:\n",
    "Plotten Sie die Dichtefunktion und Verteilungsfunktion der stetigen Gleichverteilung $U(0,1)$ mit Hilfe von `uniform` oder indem Sie die Funktionen direkt implementieren.\n",
    "\n",
    "Die Dichtefunktion der stetigen Gleichverteilung ist gegeben durch\n",
    "\n",
    "$f(x\\mid a ,b)={\\begin{cases}{\\frac  {1}{b-a}}&a\\leq x \\leq b\\\\0& \\text{sonst.}\\end{cases}}$\n",
    "\n",
    "Die zugehörige Verteilungsfunktion ist\n",
    "\n",
    "$F(x\\mid a ,b)={\\begin{cases}0&x\\leq a\\\\{\\frac  {x-a}{b-a}}&a<x<b\\\\1&x\\geq b.\\end{cases}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc28f402",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import uniform\n",
    "\n",
    "# TODO\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba418131",
   "metadata": {},
   "source": [
    "Nun betrachten wir eine gemeinsame Verteilung von zwei Zufallsvariablen $X_1,X_2$ am Beispiel der bivariaten Normalverteilung, gegeben durch\n",
    "\n",
    "${\\displaystyle f(\\mathbf {x} \\mid \\mathbf {\\mu} ,\\mathbf {\\sigma})={\\frac {1}{\\sqrt {(2\\pi )^{p}\\det(\\mathbf {\\Sigma } )}}}\\exp \\left(-{\\frac {1}{2}}({\\mathbf {x} }-{\\boldsymbol {\\mu }})^{\\top }{\\mathbf {\\Sigma } }^{-1}({\\mathbf {x} }-{\\boldsymbol {\\mu }})\\right)}$.\n",
    "\n",
    "Diese kann mit Hilfe von `multivariate_normal` geplottet werden."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6995a8d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotte bivariate Normalverteilung\n",
    "from scipy.stats import multivariate_normal\n",
    "from mpl_toolkits import mplot3d\n",
    " \n",
    "# Kovarianz setzen\n",
    "cov_val = 0.8\n",
    " \n",
    "# Erwartungswert (0,0)\n",
    "mean = np.zeros(2)\n",
    " \n",
    "# Kovarianzmatrix mit Varianz 1\n",
    "cov = np.array([[1, cov_val], [cov_val, 1]])\n",
    "     \n",
    "# Gaussverteilung mit Mean, Kovarianzmatrix\n",
    "distr = multivariate_normal(cov = cov, mean = mean)\n",
    "     \n",
    "# Plotte zwischen -3,3 mit Aufloesung 1000   \n",
    "x = np.linspace(-3, 3, 100)\n",
    "y = np.linspace(-3, 3, 100)\n",
    "X, Y = np.meshgrid(x,y)\n",
    "     \n",
    "# Speichere die Dichtefunktion und die Verteilungsfunktion\n",
    "pdf = np.zeros(X.shape)\n",
    "cdf = np.zeros(X.shape)\n",
    "for i in range(X.shape[0]):\n",
    "    for j in range(X.shape[1]):\n",
    "        pdf[i,j] = distr.pdf([X[i,j], Y[i,j]])\n",
    "        cdf[i,j] = distr.cdf([X[i,j], Y[i,j]])\n",
    "     \n",
    "# Plot Dichtefunktion\n",
    "ax = plt.axes(projection = '3d')\n",
    "ax.plot_surface(X,Y, pdf, rstride=3, cstride=3, linewidth=1, #antialiased=True,\n",
    "                cmap='viridis')\n",
    "plt.xlabel('x$_1$', fontsize=12)\n",
    "plt.ylabel('x$_2$', fontsize=12)\n",
    "ax.view_init(30, 30) #So kann man die Funktion noch drehen\n",
    "plt.show()\n",
    "\n",
    "# Plot Verteilungsfunktion\n",
    "ax = plt.axes(projection = '3d')\n",
    "ax.plot_surface(X,Y,cdf,rstride=1, cstride=1,\n",
    "                cmap='viridis', edgecolor='none')\n",
    "plt.xlabel('x$_1$', fontsize=12)\n",
    "plt.ylabel('x$_2$', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "110e7326",
   "metadata": {},
   "source": [
    "Eine Darstellung in zwei Dimensionen erfolgt mit Hilfe von Niveaumengen. Auch werden oft die Randverteilungen mit dargestellt.\n",
    "\n",
    "### Aufgabe 1.3\n",
    "Aus der zweidimensionalen Verteilung lassen sich die Randverteilungen bestimmen. Was sind die zugehoerigen Randverteilungen? Implementieren Sie dies in den vorgefertigten Funktionen `rand_f1` und `rand_f2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afea66de",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rand_f1(x):\n",
    "    # TODO\n",
    "    return x\n",
    "\n",
    "def rand_f2(x):\n",
    "    # TODO\n",
    "    return x\n",
    "\n",
    "def pdf_marg(x, y, ax, ax_histx, ax_histy):\n",
    "    ax_histx.tick_params(axis=\"x\", labelbottom=False, labeltop=True)\n",
    "    ax_histy.tick_params(axis=\"y\", labelleft=False, labelright=True)\n",
    "    ax_histx.set_ylabel('F$_1$', fontsize=12)\n",
    "    ax_histy.set_xlabel('F$_2$', fontsize=12)\n",
    "    \n",
    "    ax_histx.set_yticks([0, 0.3])\n",
    "    ax_histy.set_xticks([0, 0.3])\n",
    "    \n",
    "    # Plotte Contour\n",
    "    ax.contourf(x, y, pdf, cmap='viridis')\n",
    "\n",
    "    # Plotte Randverteilungen\n",
    "    ax_histx.plot(y, rand_f1(y),c='black', linewidth=0.5)\n",
    "    ax_histy.plot(rand_f2(y),y,c='black', linewidth=0.5)\n",
    "\n",
    "# Plot Allgemein (in einem Figure)\n",
    "fig = plt.figure()#figsize=(3, 3))\n",
    "grid = fig.add_gridspec(2, 2,  width_ratios=(4, 1), height_ratios=(1, 4),\n",
    "                      left=0.1, right=0.9, bottom=0.1, top=0.9,\n",
    "                      wspace=0.1, hspace=0.1)\n",
    "\n",
    "# Hauptachsen, die sich die Plots teilen\n",
    "ax = fig.add_subplot(grid[1, 0])\n",
    "ax_histx = fig.add_subplot(grid[0, 0], sharex=ax)\n",
    "ax_histy = fig.add_subplot(grid[1, 1], sharey=ax)\n",
    "ax.set_xlabel('x$_1$', fontsize=12)\n",
    "ax.set_ylabel('x$_2$', fontsize=12)\n",
    "\n",
    "# Plot\n",
    "pdf_marg(X,Y, ax, ax_histx, ax_histy)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0812271",
   "metadata": {},
   "source": [
    "Um die Verteilung darzustellen können wir auch Stichproben verwenden. Diese werden aus der Verteilung generiert. \n",
    "\n",
    "### Aufgabe 1.4\n",
    "Plotten Sie die Stichproben der Verteilungsfunktion. Welche Struktur ist erkennbar?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3540f0ac",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Kovarianz setzen\n",
    "cov_val = 0.8\n",
    "\n",
    "# Erwartungswert (0,0)\n",
    "mean = np.zeros(2)\n",
    "\n",
    "# Kovarianzmatrix mit Varianz 1\n",
    "cov = np.array([[1, cov_val], [cov_val, 1]]) \n",
    "\n",
    "# Gaussverteilung mit Mean, Kovarianzmatrix\n",
    "distr = multivariate_normal(cov = cov, mean = mean)\n",
    "\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeb70b4c",
   "metadata": {},
   "source": [
    "## 2 Stochastische Abhängigkeit\n",
    "\n",
    "### Aufgabe 2.1\n",
    "Schauen Sie sich die Verteilung für verschiedene Werte der Kovarianz `cov_val` zwischen -0.99 und 0.99 an. Welche Spezialfälle ergeben sich? Welche Form von Abhängigkeit sehen Sie? Wie ändern sich die Randverteilungen? Kann man aus den Randverteilungen auch die gemeinsame Verteilung bestimmen? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ed3c258",
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "816289e3",
   "metadata": {},
   "source": [
    "## 3 Copula\n",
    "\n",
    "Der Satz von Sklar zeigt auf, wie sich eine Copula mit aus der gemeinsamen Verteilungsfunktion und den Randverteilungen herleiten lässt. \n",
    "\n",
    "### Aufgabe 3.1\n",
    "Bestimmen Sie die Stichproben der Copula, die sich aus den Stichproben der obigen bivariaten Normalverteilung und den jeweiligen Randverteilungen ergibt. Benutzen Sie die Transformation auf stetig gleichverteilte Zufallsvariablen in $[0,1]$. Welcher Copula, die Sie in der Vorlesung gesehen haben, entspricht diese Verteilung?\n",
    "Wie sieht die Copula für Randverteilungen mit $N(1.0,0.5)$ aus?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3118aee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO\n",
    "#In u1,u2 sollen am Ende die Stichproben der Copula gespeichert sein\n",
    "u1=uniform.rvs(0,1,2000)\n",
    "u2=uniform.rvs(0,1,2000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec121621",
   "metadata": {},
   "source": [
    "Andersherum lässt sich mit Hilfe der Copula als Abhängigkeitsstruktur und den Randverteilungen eine gemeinsame Verteilung bestimmen. Dazu braucht man lediglich eine Transformation durch die Quantilfunktionen der Randverteilungen. Im Folgenden wird dies basierend auf der Copula aus der letzten Aufgabe und invertierten $\\chi^2$ Verteilungen mit Parameter $1$ als Randverteilungen durchgeführt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba94721f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#(1) Modelliere Copula mit Hilfe von Randverteilungen: invertierte Chi^2-Verteilung\n",
    "from scipy.stats import chi2\n",
    "#Invertierte Verteilungen sind gerade Chi^2-Verteilungen\n",
    "x1=chi2.cdf(u1,1)\n",
    "x2=chi2.cdf(u2,1)\n",
    "\n",
    "#(2)Darstellung der Verteilung mit Stichproben:\n",
    "plt.plot(x1,x2, '.', c='black', markeredgewidth = 0.1)\n",
    "plt.xlim([0.0, 0.7])\n",
    "plt.ylim([0.0, 0.7])\n",
    "plt.xlabel('x$_1$', fontsize=12)\n",
    "plt.ylabel('x$_2$', fontsize=12)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04960375",
   "metadata": {},
   "source": [
    "### Aufgabe 3.2\n",
    "Führen Sie diese Schritte für die Exponentialverteilung\n",
    "$f(x\\mid \\lambda)={\\begin{cases}{\\lambda \\exp(-\\lambda x)}& x \\geq 0 \\\\0& x<0\\end{cases}}$.\n",
    "als Randverteilungen durch. Verwenden Sie hierfür die $\\lambda=2$\n",
    "Was passiert, wenn Sie als Randverteilungen Normalverteilungen $N(0,1)$ annehmen? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "499800af",
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b7ab807",
   "metadata": {},
   "source": [
    "### Aufgabe 3.3\n",
    "Untersuchen Sie, wie sich der Korrelationskoeffizient unter den obigen Transformation verhält."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2f73d68",
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO (np.corrcoef)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6a4392e",
   "metadata": {},
   "source": [
    "## 4 Ausblick\n",
    "\n",
    "In der nächsten Vorlesung lernen wir die Konstruktion von archimedischen Copulas wie der Gumbel-Copula. Im Folgenden konstruieren wir Verteilungen aus verschiedenen archimedischen Copulas mit Randverteilungen $N(0,1)$.\n",
    "\n",
    "### 1 Gumbel-Copula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d145ff13",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Konstruiere Gumbel-Copula mit Korrelation 0.8\n",
    "def gumbel_inv_generator(u,lamb):\n",
    "    return np.exp(-u**(1./lamb))\n",
    "\n",
    "d1=np.random.uniform(0,1,3000)\n",
    "d2=np.random.uniform(0,1,3000)\n",
    "\n",
    "lamb=2.5\n",
    "\n",
    "#Sample der Verteilung mit Dichte (1-1/lamb+s/lamb)exp(-s)\n",
    "w1=np.random.gamma(1,1,3000)\n",
    "w2=np.random.gamma(1,2,3000)\n",
    "s=np.zeros(3000)\n",
    "for i in range(0,3000):\n",
    "    if(d2[i]<=lamb):\n",
    "        s[i]=w1[i]\n",
    "    else:\n",
    "        s[i]=w2[i]\n",
    "z1=d1*s**lamb\n",
    "z2=(1-d1)*s**lamb\n",
    "\n",
    "g1=gumbel_inv_generator(z1,lamb)\n",
    "g2=gumbel_inv_generator(z2,lamb)\n",
    "\n",
    "#Darstellung der Verteilung mit Stichproben:\n",
    "plt.plot(g1, g2, '.', c='black', markeredgewidth = 0.1)\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.xlabel('u$_1$', fontsize=12)\n",
    "plt.ylabel('u$_2$', fontsize=12)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect('equal', adjustable='box')\n",
    "plt.show()\n",
    "\n",
    "print(np.corrcoef(g1,g2))\n",
    "\n",
    "#Add Marginals N(0,1)\n",
    "gumb1=norm.ppf(g1)\n",
    "gumb2=norm.ppf(g2)\n",
    "\n",
    "#Darstellung der Verteilung mit Stichproben:\n",
    "plt.plot(gumb1, gumb2, '.', c='black', markeredgewidth = 0.1)\n",
    "plt.xlabel('x$_1$', fontsize=12)\n",
    "plt.ylabel('x$_2$', fontsize=12)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect('equal', adjustable='box')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eb6eee8",
   "metadata": {},
   "source": [
    "### Clayton-Copula"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e32a9cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sampling the Clayton copula\n",
    "alpha=3.2\n",
    "\n",
    "def clayton_inv_generator(u,alpha):\n",
    "    return (1+u)**(-1./alpha)\n",
    "\n",
    "#Sample from Laplace-Stieltjes-Transform\n",
    "s=np.random.gamma(1./alpha,1,3000)\n",
    "\n",
    "#transform\n",
    "c1=clayton_inv_generator(-np.log(d1)/s,alpha)\n",
    "c2=clayton_inv_generator(-np.log(d2)/s,alpha)\n",
    "\n",
    "#Darstellung der Verteilung mit Stichproben:\n",
    "plt.plot(c1, c2, '.', c='black', markeredgewidth = 0.1)\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.xlabel('u$_1$', fontsize=12)\n",
    "plt.ylabel('u$_2$', fontsize=12)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect('equal', adjustable='box')\n",
    "plt.show()\n",
    "\n",
    "print(np.corrcoef(c1,c2))\n",
    "\n",
    "#Add Marginals N(0,1)\n",
    "clay1=norm.ppf(c1)\n",
    "clay2=norm.ppf(c2)\n",
    "\n",
    "#Darstellung der Verteilung mit Stichproben:\n",
    "plt.plot(clay1, clay2, '.', c='black', markeredgewidth = 0.1)\n",
    "plt.xlabel('x$_1$', fontsize=12)\n",
    "plt.ylabel('x$_2$', fontsize=12)\n",
    "ax = plt.gca()\n",
    "ax.set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1cc5028",
   "metadata": {},
   "source": [
    "### Ausblick\n",
    "Was passiert in den Randfällen für kleine oder große Werte von $x_1$ und $x_2$? Wie verhält es sich hier bei der Gauss-Kopula?\n",
    "Wir werden hierzu die Begriffe der Tailabhängigkeit/Flankenabhängigkeit kennen lernen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6219421b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}