# -*- coding: utf-8 -*-
"""Chapter 5 exercises.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1f0JKfnLSFbJNDQhhvMqj6Yoh6GO7rbI9
"""

# Section 5.2

# Exercise 5.3

import numpy as np
from scipy.stats import norm


m = 100000
xi = 0.05
sigma = 0.25
w = 0.5
rho = [0.0, 0.5, 0.8]
eps = np.zeros((m,4))
eps[:, :2] = norm.rvs(size=(m, 2))
eps[:, 2] = eps[:,1]
eps[:, 3] = eps[:,1]
eps[:,1] = rho[0]*eps[:,0] + np.sqrt(1 - rho[0]**2)*eps[:,1]
eps[:,2] = rho[1]*eps[:,0] + np.sqrt(1 - rho[1]**2)*eps[:,2]
eps[:,3] = rho[2]*eps[:,0] + np.sqrt(1 - rho[2]**2)*eps[:,3]
R = np.exp(xi + sigma*eps)-1
Ratio1 = np.std(w*R[:,0] + (1-w)*R[:,1])/np.std(R[:,0])
Ratio2 = np.std(w*R[:,0] + (1-w)*R[:,2])/np.std(R[:,0])
Ratio3 = np.std(w*R[:,0] + (1-w)*R[:,3])/np.std(R[:,0])
Ratio = [Ratio1, Ratio2, Ratio3]

ExactRatio = [np.sqrt((1+rho[i])/2.0) for i in range(len(rho))]

print(Ratio)
print(ExactRatio)

"""# Exercise 5.4"""

import pandas as pd

url = '/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt'
x = pd.read_csv(url, header=None, sep='\s+')

x.head()

meanx = np.mean(x, axis=0)
covx = np.cov(x, rowvar=0)
R = np.exp(x)-1
meanR = np.mean(R, axis=0)
covR = np.cov(R, rowvar=0)
volR = np.sqrt(np.diag(covR))
corR = np.corrcoef(R, rowvar=0)

corR

"""# Exercise 5.6"""

a1 = 0.970
a2 = 0.767
a3 = 0.746
a4 = 0.746

a = np.array([a1,a2,a3,a4])
cor = np.outer(a, a)
np.fill_diagonal(cor, 1)

cor

"""# Section 5.3

# Exercise 5.7
"""

xi = np.array([0.006107,0.00467])
r = 0.004361
sigma = np.array([0.0158,0.019])
rho = 0.9654
eg = 0.033869

Rho = np.ones((2,2))
Rho[0][1] = rho
Rho[1][0] = rho
Sigma = np.diag(sigma) @ Rho @ np.diag(sigma)
wTilde = np.linalg.solve(Sigma,xi-r)
gamma = (eg-r)/(wTilde.T @ (xi-r))
w = gamma*wTilde
vol = np.sqrt(w.T @ Sigma @ w)

print(1- sum(w), w)

vol

rho = 0.8
Rho = np.ones((2,2))
Rho[0][1] = rho
Rho[1][0] = rho
Sigma = np.diag(sigma) @ Rho @ np.diag(sigma)
wTilde = np.linalg.solve(Sigma,xi-r)
gamma = (eg-r)/(wTilde.T @ (xi-r))
w = gamma*wTilde
vol = np.sqrt(w.T @ Sigma @ w)
print(1- sum(w), w)
print(vol)

rho = 0.6
Rho = np.ones((2,2))
Rho[0][1] = rho
Rho[1][0] = rho
Sigma = np.diag(sigma) @ Rho @ np.diag(sigma)
wTilde = np.linalg.solve(Sigma,xi-r)
gamma = (eg-r)/(wTilde.T @ (xi-r))
w = gamma*wTilde
vol = np.sqrt(w.T @ Sigma @ w)
print(1- sum(w), w)
print(vol)

"""# Exercise 5.8"""

import pandas as pd
import numpy as np

url = "/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt"
x = pd.read_csv(url, header=None, sep='\s+')

R = np.exp(x)-1
r = 0.00025
eg = np.arange(0.0005,0.0015,0.0001)
xi = np.mean(R, axis=0)
Sigma = np.cov(R, rowvar=0)

wTilde = np.linalg.solve(Sigma,xi-r)
egTilde = r + wTilde.T @ (xi-r)
volTilde = np.sqrt(wTilde.T @ Sigma @ wTilde)

print(wTilde, egTilde, volTilde)

gamma = (eg[:3]-r)/(egTilde-r)
vol = gamma*volTilde

print(gamma)
print(vol)
print(1 - gamma*sum(wTilde))

"""# Exercise 5.9"""

url = "/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt"
x = pd.read_csv(url, header=None, sep='\s+')

R = np.exp(x)-1
meanR = np.mean(R, axis=0)
covR = np.cov(R, rowvar=0)
volR = np.sqrt(np.diag(covR))
corR = np.corrcoef(R, rowvar=0)
D = np.diag(volR)
Sigma = D @ corR @ D

Sigma

"""# Exercise 5.10"""

url = "/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt"
x = pd.read_csv(url, header=None, sep='\s+')

R = np.exp(x)-1
meanR = np.mean(R, axis=0)
covR = np.cov(R, rowvar=0)
volR = np.sqrt(np.diag(covR))

r = 0.00025
xi = r + volR * np.mean((meanR-r)/volR)
wTilde = np.linalg.solve(Sigma,xi-r)
egTilde = r + wTilde.T @ (xi-r)
volTilde = np.sqrt(wTilde.T @ Sigma @ wTilde)

print(wTilde)
print(egTilde)
print(volTilde)

"""# Section 5.4

# Exercise 5.15

"""

url = "/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt"
x = pd.read_csv(url, header=None, sep='\s+')

R = np.exp(x)-1
Sigma = np.cov(R, rowvar=0)
C = np.linalg.cholesky(Sigma)

print(C)

print(C @ C.T)

print(Sigma)

xi = np.mean(R, axis=0)
m = 100000

eps = np.random.normal(size=(m,4))
Rstar = (np.tile(xi, (m,1)) + eps @ C.T)

print(xi)
print(np.mean(Rstar, axis=0))

print(np.sqrt(np.diag(Sigma)))
print(np.std(Rstar, axis=0))

np.diag(np.linalg.solve(np.diag(np.sqrt(np.diag(Sigma))), np.ones(4))) @ Sigma @ np.diag(np.linalg.solve(np.diag(np.sqrt(np.diag(Sigma))), np.ones(4)))

np.corrcoef(Rstar, rowvar=0)

"""# Exercise 5.16"""

url = "/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt"
x = pd.read_csv(url, header=None, sep='\s+')

R = np.exp(x)-1
Sigma = np.cov(R, rowvar=0)
C = np.linalg.cholesky(Sigma)

xi = np.array(np.mean(R, axis=0))
xi = xi.reshape(-1,1)
w = np.ones(4) * 0.25
m = 100000

eps = np.random.normal(size=(4,m))
Rstar = (xi+C @ eps).T

Rcal = Rstar @ w

print(np.std(Rcal))
print(np.std(Rstar, axis=0))

print(np.corrcoef(Rstar, rowvar=0))

"""# Exercise 5.18"""

ER = 0.11
sdR = 0.25
corR1 = 0.8
corR2 = 0.4

sigma = np.sqrt(np.log(1+sdR**2/((1+ER)**2)))
xi = np.log(1+ER)-0.5*sigma**2
rho1 = np.log(1+sdR**2*corR1/((1+ER)**2))/sigma**2
rho2 = np.log(1+sdR**2*corR2/((1+ER)**2))/sigma**2

print(xi,sigma,rho1,rho2)

"""# Exercise 5.20"""

from scipy.stats import gamma

url = "/studier/emner/matnat/math/STK3505/h20/newyorkdaily.txt"
x = pd.read_csv(url, header=None, sep='\s+')

R = np.exp(x)-1

m = 100000
alpha = 500
Sigma = np.cov(R, rowvar=0)
C = np.linalg.cholesky(Sigma)
xi = np.array(np.mean(R, axis=0))
xi = xi.reshape(-1,1)

eta = np.random.normal(size=(4,m))
eps = (C @ eta).T
Z = (alpha - 1)*gamma.rvs(alpha, size=(m,1))
X = xi + (np.sqrt(Z)*eps).T
mean = np.mean(X, axis=1)
cov = np.cov(X, rowvar=1)

cov

"""# Exercise 5.21"""

Rcal = np.mean(np.exp(X)-1, axis=0)
ind = [0.01,0.10,0.25,0.50,0.75,0.90,0.99]
percentile = [np.sort(Rcal)[int(m*i)] for i in ind]

percentile

"""# Section 5.5

# Exercise 5.22
"""

from scipy.stats import norm

xi = 0.05
sigma = 0.25
K = 20
yps = np.array([0.05,0.25,0.5,0.75,0.95])

print(np.exp(K*xi+np.sqrt(K)*sigma*norm.ppf(yps))-1)

"""# Exercies 5.23"""

import matplotlib.pyplot as plt

xi = 0.05
sigma = 0.25
K = 20
m = 50

eta = norm.rvs(size=(K,m))
R0k = np.exp(np.cumsum(xi+sigma*eta, axis=0))-1
R0k = np.insert(R0k, 0, np.zeros(m), axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), R0k)

"""# Exercise 5.24"""

from scipy.stats import gamma

xi = 0.05
sigma = 0.25
xisigma = 0.25
alpha = 5
K = 20
m = 50

eta = norm.rvs(size=(K, m))
Gk = gamma.rvs(alpha,size=(K,m), scale=1/alpha)
sigmak = xisigma/np.sqrt(Gk)
R0k = np.exp(np.cumsum(xi+sigmak*eta,axis=0))-1
R0k = np.insert(R0k, 0, np.zeros(m), axis=0)

for i in range(m):
  plt.plot(np.arange(K+1), R0k)

"""# Exercise 5.25"""

xi = 0.05
sigma = 0.25
xisigma = 0.25
alpha = 5
K = 20
yps = [0.05,0.25,0.5,0.75,0.95]
m = 10000

eta = norm.rvs(size=(K,m))
Gk = gamma.rvs(alpha,size=(K,m), scale=1/alpha)
sigmak = xisigma/np.sqrt(Gk)
R0K = np.exp(np.sum(xi+sigmak*eta,axis=0))-1

print([np.sort(R0K)[int(y*m)] for y in yps])
print([np.exp(K*xi+np.sqrt(K)*sigma*norm.ppf(y))-1 for y in yps])

"""# Exercise 5.26"""

sigma = 0.25
K = 1000000

eta = norm.rvs(size=K)
X0k = sigma*np.cumsum(eta)
ind = np.arange(100,K,100)

plt.plot(ind,X0k[ind])

"""# Exercise 5.27"""

import numpy as np
import matplotlib.pyplot as plt

xi = 0.05
sigma = 0.25
rho1 = 0
rho2 = -0.9
rho3 = 0.9
K = 50

eta = np.random.normal(size=(K,2))

eps1 = np.copy(eta)
eps1[:,1] = rho1*eps1[:,0]+np.sqrt(1-rho1**2)*eps1[:,1]
R0k1 = np.exp(np.cumsum(xi+sigma*eps1,axis=0))-1
R0k1 = np.insert(R0k1, 0, np.zeros(2), axis=0)

eps2 = np.copy(eta)
eps2[:,1] = rho2*eps2[:,0]+np.sqrt(1-rho2**2)*eps2[:,1]
R0k2 = np.exp(np.cumsum(xi+sigma*eps2,axis=0))-1
R0k2 = np.insert(R0k2, 0, np.zeros(2), axis=0)

eps3 = np.copy(eta)
eps3[:,1] = rho3*eps3[:,0]+np.sqrt(1-rho3**2)*eps3[:,1]
R0k3 = np.exp(np.cumsum(xi+sigma*eps3,axis=0))-1
R0k3 = np.insert(R0k3, 0, np.zeros(2), axis=0)

fig, axs = plt.subplots(3)
ymin = min(np.min(R0k1),np.min(R0k2),np.min(R0k3))
ymax = max(np.max(R0k1),np.max(R0k2),np.max(R0k3))
axs[0].plot(R0k1)
axs[0].axis(ymin=ymin-20,ymax=ymax+20)
axs[1].plot(R0k2)
axs[1].axis(ymin=ymin-20,ymax=ymax+20)
axs[2].plot(R0k3)
axs[2].axis(ymin=ymin-20,ymax=ymax+20)

"""# Section 5.6

# Exercise 5.28
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics import tsaplots

url = "/studier/emner/matnat/math/STK3505/h20/nibornor.txt"
r = pd.read_csv(url, header=None, sep='\s+')

plt.scatter(np.arange(len(r)),r)
plt.xticks(np.arange(len(r)))

autocor = sm.tsa.acf(r)
a = autocor[1]

plt.scatter(np.arange(len(autocor)), autocor)
plt.plot(np.arange(len(autocor)), np.power(a, np.arange(len(autocor))))
plt.hlines(0.0, xmin=0, xmax=20)
plt.xticks(np.arange(len(r)))

a

"""# Exercise 5.29"""

url = "/studier/emner/matnat/math/STK3505/h20/nibornor.txt"
r = pd.read_csv(url, header=None, sep='\s+')

xiHat = np.mean(r)
autocor = sm.tsa.acf(r)
aHat = autocor[1]
sigmaHat = np.std(r)*np.sqrt(1-aHat**2)

print(*xiHat)
print(aHat)
print(*sigmaHat)

"""# Exercise 5.30"""

xi = 0.062
sigma = 0.018
a = 0.71
r0 = 0.04
K = 25
m = 25

x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi

plt.plot(np.arange(K+1), r)

r0 = 0.1
x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi
plt.plot(np.arange(K+1), r)

"""# Exercise 5.31"""

xi = 0.062
sigma = 0.018
a = 0.71
r0 = 0.04
K = 10
m = 100000

x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi

from scipy.stats import kde

for k in range(1, K+1):
  density = kde.gaussian_kde(r[k,:])
  x = np.linspace(-0.05, 0.15, 200)
  y = density(x)
  plt.plot(x, y)

r0 = 0.1
x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi
for k in range(1, K+1):
  density = kde.gaussian_kde(r[k,:])
  x = np.linspace(-0.05, 0.15, 200)
  y = density(x)
  plt.plot(x, y)

"""# Exercise 5.32"""

xi = 0.04
sigma = 0.018
a = 0.71
r0 = 0.04
K = 25
m = 25

x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi
plt.plot(np.arange(K+1), r)

r0 = 0.1
x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi
plt.plot(np.arange(K+1), r)
plt.hlines(0.0, xmin=0, xmax=25)

"""# Exercise 5.35"""

xi = 0.062
sigma = 0.28
a = 0.71
sigmax = sigma/np.sqrt(1-a**2)
r0 = 0.04
K = 25
m = 25

x0 = np.log(r0/xi)+0.5*sigmax**2
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = xi*np.exp(-0.5*sigmax**2+X)
plt.plot(np.arange(K+1), r)

"""# Exercise 5.36"""

xi = 0.062
sigma = 0.018
a = 0.71
r0 = 0.04
K = 20
m = 100000

x0 = r0-xi
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = X+xi
r0K = np.prod(1+r, axis=0) - 1

density = kde.gaussian_kde(r0K)
x = np.linspace(0.0, 10.0, 200)
y = density(x)
plt.plot(x, y)
plt.ylim(0.0, 1.0)

sigma = 0.18
sigmax = sigma/np.sqrt(1-a**2)

x0 = np.log(r0/xi)+0.5*sigmax**2
X = np.ones((K+1,m))*x0
for k in range(K):
  eps = np.random.normal(size=m)
  X[k+1,:] = a*X[k,:]+sigma*eps
r = xi*np.exp(-0.5*sigmax**2+X)
r0K = np.prod(1+r, axis=0) - 1

density = kde.gaussian_kde(r0K)
x = np.linspace(0.0, 10.0, 200)
y = density(x)
plt.plot(x, y)
plt.ylim(0.0, 1.0)