#!/usr/bin/env python
# coding: utf-8

# <h1> Jupyter-Notebook zu Copulas (Musterlösung) </h1>
# (Lisa Kusch, Nov 2022)

# ## 1 Grundlagen: Wahrscheinlichkeitsdichte und Wahrscheinlichkeitsverteilung
# 
# Die Dichtefunktion der Normalverteilung ist gegeben durch
# 
# $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 $.
# 
# Die zugehörige Verteilungsfunktion ist
# 
# $F(x\mid \mu, \sigma)= {\frac {1}{2}}\left(1+\operatorname {erf} \left({\frac {x-\mu}{\sqrt {2}\sigma}}\right)\right)$,
# 
# wobei $\operatorname{erf}$ die Fehlerfunktion ist.
# 
# 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.

# In[ ]:


# Benötigte Imports
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.stats import norm

# Seed wird festgesetzt (Reproduzierbarkeit)
np.random.seed(1)


# Nun wird die Funktion im Bereich $[-4,4]$ geplottet. Dabei ist $\mu = 0$ und $\sigma = 1$. 

# In[ ]:


# Plotte Normalverteilung mit Erwartungswert mu und Standardabweichung sigma
mu=0.0
sigma=1.0

#Plotte zwischen -4 und 4 mit Auflösung von 0.005 (Schrittweite)
x = np.arange(-4, 4, 0.005)

plt.plot(x, norm.pdf(x,mu,sigma), c='black')
plt.xlabel('x', fontsize=13)
plt.ylabel('f(x)', fontsize=13)
plt.show()

plt.plot(x, norm.cdf(x,mu,sigma), c='black')
plt.xlabel('x', fontsize=13)
plt.ylabel('F(x)', fontsize=13)
plt.show()
    


# ### Aufgabe 1.1:
# Probieren Sie verschiedene Werte von $\mu$ und $\sigma$ aus. Wie verändert sich die Funktion? Was bedeuten die beiden Parameter?
# 
# ### Aufgabe 1.2:
# Plotten Sie die Dichtefunktion und Verteilungsfunktion der stetigen Gleichverteilung $U(0,1)$ mit Hilfe von `uniform` oder indem Sie die Funktionen direkt implementieren.
# 
# Die Dichtefunktion der stetigen Gleichverteilung ist gegeben durch
# 
# $f(x\mid a ,b)={\begin{cases}{\frac  {1}{b-a}}&a\leq x \leq b\\0& \text{sonst.}\end{cases}}$
# 
# Die zugehörige Verteilungsfunktion ist
# 
# $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}}$

# In[ ]:


from scipy.stats import uniform

# Plotte stetige Gleichverteilung U(0,1) auf [0,1]
a=0.0
b=1.0

#Plotte zwischen -0.5 und 1.5 mit Auflösung von 0.005 (Schrittweite)
x = np.arange(-0.5, 1.5, 0.005)
x_new=np.copy(x)
x_new[100]=np.nan
x_new[300]=np.nan

plt.plot(x_new, uniform.pdf(x_new,a,b), c='black')
plt.xlabel('x', fontsize=13)
plt.ylabel('f(x)', fontsize=13)
plt.show()

plt.plot(x, uniform.cdf(x,a,b), c='black')
plt.xlabel('x', fontsize=13)
plt.ylabel('F(x)', fontsize=13)
plt.show()
    


# Nun betrachten wir eine gemeinsame Verteilung von zwei Zufallsvariablen $X_1,X_2$ am Beispiel der bivariaten Normalverteilung, gegeben durch
# 
# ${\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)}$.
# 
# Diese kann mit Hilfe von `multivariate_normal` geplottet werden.

# In[ ]:


# Plotte bivariate Normalverteilung
from scipy.stats import multivariate_normal
from mpl_toolkits import mplot3d
 
# Kovarianz setzen
cov_val = 0.8
 
# Erwartungswert (0,0)
mean = np.zeros(2)
 
# Kovarianzmatrix mit Varianz 1
cov = np.array([[1, cov_val], [cov_val, 1]])
     
# Gaussverteilung mit Mean, Kovarianzmatrix
distr = multivariate_normal(cov = cov, mean = mean)
     
# Plotte zwischen -3,3 mit Aufloesung 1000   
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
X, Y = np.meshgrid(x,y)
     
# Speichere die Dichtefunktion und die Verteilungsfunktion
pdf = np.zeros(X.shape)
cdf = np.zeros(X.shape)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        pdf[i,j] = distr.pdf([X[i,j], Y[i,j]])
        cdf[i,j] = distr.cdf([X[i,j], Y[i,j]])
     
# Plot Dichtefunktion
ax = plt.axes(projection = '3d')
ax.plot_surface(X,Y, pdf, rstride=3, cstride=3, linewidth=1, #antialiased=True,
                cmap='viridis')
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax.view_init(30, 30) #So kann man die Funktion noch drehen
plt.show()

# Plot Verteilungsfunktion
ax = plt.axes(projection = '3d')
ax.plot_surface(X,Y,cdf,rstride=1, cstride=1,
                cmap='viridis', edgecolor='none')
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
plt.show()


# Eine Darstellung in zwei Dimensionen erfolgt mit Hilfe von Niveaumengen. Auch werden oft die Randverteilungen mit dargestellt.
# 
# ### Aufgabe 1.3
# 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`.

# In[ ]:


def rand_f1(x):
    return norm.pdf(x,0,1)

def rand_f2(x):
    return norm.pdf(x,0,1)

def pdf_marg(x, y, ax, ax_histx, ax_histy):
    ax_histx.tick_params(axis="x", labelbottom=False, labeltop=True)
    ax_histy.tick_params(axis="y", labelleft=False, labelright=True)
    ax_histx.set_ylabel('F$_1$', fontsize=12)
    ax_histy.set_xlabel('F$_2$', fontsize=12)
    
    ax_histx.set_yticks([0, 0.3])
    ax_histy.set_xticks([0, 0.3])
    
    # Plotte Contour
    ax.contourf(x, y, pdf, cmap='viridis')

    # Plotte Randverteilungen
    ax_histx.plot(y, rand_f1(y),c='black', linewidth=0.5)
    ax_histy.plot(rand_f2(y),y,c='black', linewidth=0.5)

# Plot Allgemein (in einem Figure)
fig = plt.figure()#figsize=(3, 3))
grid = fig.add_gridspec(2, 2,  width_ratios=(4, 1), height_ratios=(1, 4),
                      left=0.1, right=0.9, bottom=0.1, top=0.9,
                      wspace=0.1, hspace=0.1)

# Hauptachsen, die sich die Plots teilen
ax = fig.add_subplot(grid[1, 0])
ax_histx = fig.add_subplot(grid[0, 0], sharex=ax)
ax_histy = fig.add_subplot(grid[1, 1], sharey=ax)
ax.set_xlabel('x$_1$', fontsize=12)
ax.set_ylabel('x$_2$', fontsize=12)

# Plot
pdf_marg(X,Y, ax, ax_histx, ax_histy)
plt.show()


# Um die Verteilung darzustellen können wir auch Stichproben verwenden. Diese werden aus der Verteilung generiert. 
# 
# ### Aufgabe 1.4
# Plotten Sie die Stichproben der Verteilungsfunktion. Welche Struktur ist erkennbar?

# In[ ]:


# Kovarianz setzen
cov_val = 0.8
 
# Erwartungswert (0,0)
mean = np.zeros(2)
 
# Kovarianzmatrix mit Varianz 1
cov = np.array([[1, cov_val], [cov_val, 1]])
     
# Gaussverteilung mit Mean, Kovarianzmatrix
distr = multivariate_normal(cov = cov, mean = mean)

# Generiere 2000 Stichproben
data = distr.rvs(size = 2000)
     
# Plot
plt.plot(data[:,0],data[:,1], '.', c='black', markeredgewidth = 0.1)
plt.xlim([-4, 4])
plt.ylim([-4, 4])
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box') #Setze Achsen gleichlang
plt.show()

#Wir veraendern die Kovarianz cov_val: 0: Kreis mit gleichem Radius, 1: perfekte lineare Abhaengikeit, -1: perfekte lineare Abh. im negativen Sinne
#Randverteilung bleibt gleich


# ## 2 Stochastische Abhängigkeit
# 
# ### Aufgabe 2.1
# 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? 

# In[ ]:


#Wir veraendern die Kovarianz cov_val: 0: Kreis mit gleichem Radius, 1: perfekte lineare Abhaengikeit, -1: perfekte lineare Abh. im negativen Sinne
#Randverteilung bleibt gleich


# ## 3 Copula
# 
# Der Satz von Sklar zeigt auf, wie sich eine Copula mit aus der gemeinsamen Verteilungsfunktion und den Randverteilungen herleiten lässt. 
# 
# ### Aufgabe 3.1
# 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?
# Wie sieht die Copula für Randverteilungen mit $N(1.0,0.5)$ aus?

# In[ ]:


#(1) Transformation mittels Randverteilungen F_1, F_2 (normalverteilt mit N(0,1))
u1=norm.cdf(data[:,0],0,1)
u2=norm.cdf(data[:,1],0,1)

#(2) Darstellung der Copula mit Stichproben:
plt.plot(u1,u2, '.', c='black', markeredgewidth = 0.1)
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.xlabel('u$_1$', fontsize=12)
plt.ylabel('u$_1$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()

#Gauss-Copula


# 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.

# In[ ]:


#(1) Modelliere Copula mit Hilfe von Randverteilungen: invertierte Chi^2-Verteilung
from scipy.stats import chi2
#Invertierte Verteilungen sind gerade Chi^2-Verteilungen
x1=chi2.cdf(u1,1)
x2=chi2.cdf(u2,1)

#(2)Darstellung der Verteilung mit Stichproben:
plt.plot(x1,x2, '.', c='black', markeredgewidth = 0.1)
plt.xlim([0.0, 0.7])
plt.ylim([0.0, 0.7])
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()


# ### Aufgabe 3.2
# Führen Sie diese Schritte für die Exponentialverteilung
# $f(x\mid \lambda)={\begin{cases}{\lambda \exp(-\lambda x)}& x \geq 0 \\0& x<0\end{cases}}$.
# als Randverteilungen durch. Verwenden Sie hierfür die $\lambda=2$
# Was passiert, wenn Sie als Randverteilungen Normalverteilungen $N(0,1)$ annehmen? 

# In[ ]:


#(1) Modelliere Copula mit Hilfe von Randverteilungen: Exponentialverteilung
def inv_exp(x,mu):
    return -np.log(1-x)/mu

y1=inv_exp(u1,2)
y2=inv_exp(u2,2)

#(2) Darstellung der Verteilung mit Stichproben:
plt.plot(y1,y2, '.', c='black', markeredgewidth = 0.1)
plt.xlim([-0.1, 4])
plt.ylim([-0.1, 4])
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()

# Modelliere Copula mit Hilfe von Randverteilungen: Normalverteilung
from scipy.stats import invgauss

z1=norm.ppf(u1)
z2=norm.ppf(u2)

#Darstellung der Verteilung mit Stichproben:
plt.plot(z1,z2, '.', c='black', markeredgewidth = 0.1)
plt.xlim([-4, 4])
plt.ylim([-4, 4])
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()


# ### Aufgabe 3.3
# Untersuchen Sie, wie sich der Korrelationskoeffizient unter den obigen Transformation verhält.

# In[ ]:


#Korrelation ändert sich mit Transformation
print(np.corrcoef(data[:,0], data[:,1]))
print(np.corrcoef(u1, u2))
print(np.corrcoef(x1, x2))
print(np.corrcoef(y1, y2))


# ## 4 Ausblick
# 
# 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)$.
# 
# ### 1 Gumbel-Copula

# In[ ]:


#Konstruiere Gumbel-Copula mit Korrelation 0.8
def gumbel_inv_generator(u,lamb):
    return np.exp(-u**(1./lamb))

d1=np.random.uniform(0,1,3000)
d2=np.random.uniform(0,1,3000)

lamb=2.5

#Sample der Verteilung mit Dichte (1-1/lamb+s/lamb)exp(-s)
w1=np.random.gamma(1,1,3000)
w2=np.random.gamma(1,2,3000)
s=np.zeros(3000)
for i in range(0,3000):
    if(d2[i]<=lamb):
        s[i]=w1[i]
    else:
        s[i]=w2[i]
z1=d1*s**lamb
z2=(1-d1)*s**lamb

g1=gumbel_inv_generator(z1,lamb)
g2=gumbel_inv_generator(z2,lamb)

#Darstellung der Verteilung mit Stichproben:
plt.plot(g1, g2, '.', c='black', markeredgewidth = 0.1)
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.xlabel('u$_1$', fontsize=12)
plt.ylabel('u$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()

print(np.corrcoef(g1,g2))

#Add Marginals N(0,1)
gumb1=norm.ppf(g1)
gumb2=norm.ppf(g2)

#Darstellung der Verteilung mit Stichproben:
plt.plot(gumb1, gumb2, '.', c='black', markeredgewidth = 0.1)
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()


# ### Clayton-Copula

# In[ ]:


#Sampling the Clayton copula
alpha=3.2

def clayton_inv_generator(u,alpha):
    return (1+u)**(-1./alpha)

#Sample from Laplace-Stieltjes-Transform
s=np.random.gamma(1./alpha,1,3000)

#transform
c1=clayton_inv_generator(-np.log(d1)/s,alpha)
c2=clayton_inv_generator(-np.log(d2)/s,alpha)

#Darstellung der Verteilung mit Stichproben:
plt.plot(c1, c2, '.', c='black', markeredgewidth = 0.1)
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.xlabel('u$_1$', fontsize=12)
plt.ylabel('u$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()

print(np.corrcoef(c1,c2))

#Add Marginals N(0,1)
clay1=norm.ppf(c1)
clay2=norm.ppf(c2)

#Darstellung der Verteilung mit Stichproben:
plt.plot(clay1, clay2, '.', c='black', markeredgewidth = 0.1)
plt.xlabel('x$_1$', fontsize=12)
plt.ylabel('x$_2$', fontsize=12)
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.show()


# ### Ausblick
# 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?
# Wir werden hierzu die Begriffe der Tailabhängigkeit/Flankenabhängigkeit kennen lernen.