0
点赞
收藏
分享

微信扫一扫

Codes for VaR-Copula

Silence潇湘夜雨 2022-01-31 阅读 20

Codes (Latex Environment)

\chapter{附件}

\section{Python代码}

编程软件:Python3.8\ jupyter\ notebook

电脑系统:Window\ 10\ 家庭中文版

电脑型号:ASUS\ TUF\ Gaming\ FA506IU

处理器:AMD\ Ryzen\ 7\ 4800H\ with\ Radeon\ Graphics\ 2.90GHz

电脑内存(RAM):16.0GB

\subsection{收益率边缘分布}

\begin{lstlisting}[language=Python]

import numpy as np

import pandas as pd

import scipy.stats as st

from scipy.stats import t

import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

sh = pd.read_csv('上证指数历史数据(2).csv')

sz = pd.read_csv('深证成份指数历史数据(2).csv')

hs = pd.read_csv('香港恒生指数历史数据.csv')

def get_pct(data):

    data.index = data['日期']

    data_pct = data ['涨跌幅']

    return data_pct

sh_pct = get_pct(sh)

sz_pct = get_pct(sz)

hs_pct = get_pct(hs)

#边缘密度

def show_edge(data,choose):

    mu = data.mean()#.values

    sigma = data.std()#.values

    num_bins = 100

    x = mu + sigma * np.random.randn(1000000)

    # 正态分布的数据

    data.hist(bins=300,figsize=(16,8),density = True,range=(-0.1,0.1)

                 ,color='b',alpha=0.6,label='密度函数')

    n, bins, patches = plt.hist(x, num_bins, density = True,alpha=0)

    # 拟合曲线

    y = st.norm.pdf(bins, mu, sigma)

    plt.plot(bins, y,'r', alpha = 0.8,linewidth=2.5,label='正态分布')

    #data.plot(kind = 'kde', color = 'g',linewidth=1.5,alpha=1,

                  label = '核密度图')

    t_dimen = 10

    z = np.linspace(t.ppf(0.01, t_dimen,loc=mu,scale=sigma),t.ppf(0.99

                 ,t_dimen,loc=mu,scale=sigma), 1000)

    plt.plot(z, t.pdf(z, t_dimen,loc=mu,scale=sigma), alpha=1.0

                 ,linewidth=2.5, c='gray',label='t-分布')

    scipy_kde=st.gaussian_kde(data)#高斯核密度估计

    X_plot=np.linspace(-0.1,0.1,1000)

    dens=scipy_kde.evaluate(X_plot)

    plt.plot(X_plot,dens,c='g',linewidth=2.5,label='核密度图')

    font = {'size': 15}#设置横纵坐标的名称以及对应字体格式、大小

    plt.xlabel('Expectation',font)

    plt.ylabel('Probability',font)

    plt.title('histogram of normal distribution: $\mu =$' +

                 str(mu)[:6] + ', $\sigma =$'+ str(sigma)[:6],font)

    axes = plt.gca()

    axes.set_xlim([-0.075,0.075])

    #显示题注

    plt.legend(loc=0)

    plt.legend(fontsize=15)#显示图例,设置图例字体大小

    plt.grid(True)

    plt.savefig(str(choose)+'.png')

    plt.show()

show_edge(sh_pct,choose='sh')

show_edge(sz_pct,choose='sz')

show_edge(hs_pct,choose='hs')

time_preiod = ['2015年12月9日','2016年12月9日','2017年12月8日'

        ,'2018年12月7日','2019年12月9日','2020年12月9日']

def get_period_pct(time_preiod):

    n = len(time_preiod)

    S=[]

    for i in range(n-1):

        sh_mean,sh_std = sh_pct[time_preiod[i]:time_preiod[i+1]]

                 .mean(),sh_pct[time_preiod[i]:time_preiod[i+1]].std()

        sz_mean,sz_std = sz_pct[time_preiod[i]:time_preiod[i+1]]

                 .mean(),sz_pct[time_preiod[i]:time_preiod[i+1]].std()

        hs_mean,hs_std = hs_pct[time_preiod[i]:time_preiod[i+1]]

                 .mean(),hs_pct[time_preiod[i]:time_preiod[i+1]].std()

        S.append([sh_mean,sh_std,sz_mean,sz_std,hs_mean,hs_std])

    return pd.DataFrame(S)

get_period_pct(time_preiod)

\end{lstlisting}

\subsection{相关系数及等高线}

\begin{lstlisting}[language=Python]

def merge(datalist):

    n=len(datalist)

    result = pd.DataFrame(datalist[0]).copy()

    result.columns=[0]

    for i in range(1,n):

        result[i]=datalist[i]

    result = result.dropna()

    return result

a= [sh_pct,sz_pct,hs_pct]

data = merge(a)

data.corr(method='pearson', min_periods=1)

data.corr(method='kendall', min_periods=1)

data.corr(method='spearman', min_periods=1)

#分布图

def show_mult(data,a,b,choose):

    plt.figure(figsize=(16,8))

    font = {'size': 15}#设置横纵坐标的名称以及对应字体格式、大小

    plt.scatter(data.values[:,a],data.values[:,b],marker = '+')

    plt.xlabel('深证指数',font)

    plt.ylabel('恒生指数',font)

    plt.title('scatter figure',font)

    plt.grid(True)

    plt.savefig(str(choose)+'.png')

    plt.show()

#这里修改参数即可得到不同数据集之间散点图,下同

show_mult(data,1,2,choose='scatter3')

#等高线

def norm_d(x1,x2,mu,Sigma):

    mu1 = mu[0]

    mu2 = mu[1]

    s1,s2 = Sigma[0,0]**0.5,Sigma[1,1]**0.5

    rho = Sigma[0,1]/s1/s2

    num=((1)/(2*np.pi*s1*s2*np.sqrt(1-rho**2)))

    A=((x1-mu1)**2)/(s1**2)

    B=2*rho*(((x1-mu1)*(x2-mu2))/(s1*s2))

    C=((x2-mu2)**2)/(s2**2)

    D=-1/(2*(1-rho**2))*(A-B+C)

    pdf=num*np.exp(D)

    return pdf

def show_contour(data,a,b,choose):

    data1 = data[[a,b]]

    plt.style.use('seaborn-whitegrid')

    x1=np.arange(-0.1,0.1,0.001)

    y1=np.arange(-0.1,0.1,0.001)

    x1,y1 = np.meshgrid(x1,y1)

    plt.figure(figsize=(16,8))

    mu = data1.mean().values

    Sigma = np.sqrt(data1.cov().values)

    z = norm_d(x1,y1,mu,Sigma)

    plt.contourf(x1, y1, z, 10)

    plt.scatter(data.values[:,a],data.values[:,b],marker = 'x')

    plt.grid(True)

    font = {'size': 15}#设置横纵坐标的名称以及对应字体格式、大小

    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']

    plt.xlabel('上证指数',font)

    plt.ylabel('深证指数',font)

    plt.title('contour figure',font)

    plt.savefig(str(choose)+'.png')

    plt.show()

show_contour(data,0,1,choose='contour1')

\end{lstlisting}

\subsection{VaR计算}

\begin{lstlisting}[language=Python]

import scipy.stats as st

def get_VaR(data,choose):

    if choose == 0:

        return np.ptp(data),st.skew(data),st.kurtosis(data)

                         ,np.percentile(data,[5,10,90,95])

    if choose ==1:

        mu = data.mean()

        sigma = data.std()

        num_bins = 100

        x = mu + sigma * np.random.randn(1000000)

        return np.percentile(x,[5,10,90,95])

    if choose ==2:

        mu = data.mean()

        sigma = data.std()

        X = st.t(df=10,loc=mu, scale=sigma)

        return X.ppf(0.05),X.ppf(0.1),X.ppf(0.9),X.ppf(0.95)

for i in range(3):

        get_VaR(sh_pct,i)

        get_VaR(sz_pct,i)

        get_VaR(hs_pct,i)

def MCMC_kde(data,choose):

    scipy_kde=st.gaussian_kde(data)

    X_plot=np.linspace(-0.12,0.12,1201)

    dens=scipy_kde.evaluate(X_plot)

    X=[]

    x=0

    for i in range(100000):

        z1 = np.random.randint(-100,100)/5000

        y=x+z1

        pro = dens[int((y+0.12)*5000)]/dens[int((x+0.12)*5000)]

        z2 = np.random.uniform()

        if pro >z2:

            X.append(y)

            x=y

        else:

            X.append(x)

    data_MCMC = pd.DataFrame(X[50000:100000])

    data_MCMC.hist(bins=300,figsize=(16,8),density = True

                 ,range=(-0.1,0.1),color='b',alpha=0.6,label='MCMC抽样')

    plt.plot(X_plot,dens,c='r',linewidth=2.5,label='核密度值')

    font = {'size': 15}#设置横纵坐标的名称以及对应字体格式、大小

    plt.xlabel('Expectation',font)

    plt.ylabel('Probability',font)

    plt.title('MCMC取样',font)

    axes = plt.gca()

    axes.set_xlim([-0.075,0.075])

    #显示题注

    plt.legend(loc=0)

    plt.legend(fontsize=15)#显示图例,设置图例字体大小

    plt.grid(True)

    plt.savefig(str(choose)+'MCMC.png')

    plt.show()

    print(np.percentile(data_MCMC,[5,10,90,95]))

MCMC_kde(sh_pct,choose='sh_pct')

MCMC_kde(sz_pct,choose='sz_pct')

MCMC_kde(hs_pct,choose='hs_pct')

\end{lstlisting}

\subsection{相依性计算}

\begin{lstlisting}[language=Python]

def rank(df):

    df = df.rank(method='min',axis=0)

    df=df.T

    return (df)

def pct_to_rank(data0):

    data=data0.copy()

    col = data.columns

    n = len(col)

    for i in range(n):

        data[col[i]]=rank(data[col[i]])

    return data

data_rank = pct_to_rank(data)

def depends(data,a,precent,b):

    df = data[[a,b]].copy()

    n = len(df.index)

    k=0

    for i in range(n):

        if df[a][i]>=n*precent and df[b][i]>=n*precent and precent>=0.5:

            k = k+1

        if df[a][i]<=n*precent and df[b][i]<=n*precent and precent<0.5:

            k = k+1

    return k/(n*precent)

def show_depends(data,a,b,start,end,n):

    l = []

    step = (end-start)/n

    for i in range(n+1):

        precent = start+i*step

        l.append(depends(data,a,precent,b))

    plt.figure(figsize=(16,8))

    plt.plot(l)

    plt.grid(True)

    plt.show()

show_depends(data_rank,0,1,0.9,1,10)

#注意到a不能取0

show_depends(data_rank,0,1,0.01,0.11,10)

\end{lstlisting}

\subsection{MCMC模拟}

\begin{lstlisting}[language=Python]

import scipy

import sympy

from sympy import diff, symbols

import numpy as np

import pandas as pd

import math

from math import *

import matplotlib.pyplot as plt

x = symbols('x', real=True)

y = symbols('y', real=True)

z = symbols('z', real=True)

e = sympy.exp(1)

t = symbols('t', real=True)

def fun_diff3(function,x,y,z):

function1 = diff(function, x, 1)

function2 = diff(function1, y, 1)

function3 = diff(function2, z, 1)

return (function3)

def pdf(sth_pdf,x0,y0,z0,t0):

if 0<x0<1 and 0<y0<1 and 0<z0<1:

return float(sth_pdf.subs({x:x0,y:y0,z:z0,t:t0,e:math.e}))

else:

return 0

def get_MCMC_values(func_pdf,t0,itn):

L=[]

l=np.array([0.5,0.5,0.5])

for i in range(itn):

z1 = (np.random.random()-0.5)/5

z2 = (np.random.random()-0.5)/5

z3 = (np.random.random()-0.5)/5

u = np.random.random()

z4 = pdf(func_pdf,l[0],l[1],l[2],t0)

l_new = l+np.array([z1,z2,z3])

z5 = pdf(func_pdf,l_new[0],l_new[1],l_new[2],t0)

if (z5/z4)>u:

L.append(l_new)

l=l_new

else:

L.append(l)

return (pd.DataFrame(L))

def get_VaR(L,precent):

S=[]

n=len(L)

for i in range(n):

S.append(L.iloc[i].sum())

S=np.array(S)

return np.percentile(S,precent)

gumble_cdf = e**(-((-sympy.log(x))**t+(-sympy.log(y))**t+(-sympy.log(z))**t)**(1/t))

gumble_pdf = fun_diff3(gumble_cdf,x,y,z)

gumble_values = get_MCMC_values(gumble_pdf,t0=1.788,itn=2000)

print(get_VaR(gumble_values,[5,10,90,95]))

clayton_cdf = (x**(-t)+y**(-t)+z**(-t)-2)**(-1/t)

clayton_pdf = fun_diff3(clayton_cdf,x,y,z)

clayton_values = get_MCMC_values(clayton_pdf,t0=1.136,itn=2000)

print(get_VaR(clayton_values,[5,10,90,95]))

\end{lstlisting}

\section{R代码}

\subsection{copula拟合}

\begin{lstlisting}[]

library(parallel)

library(zoo)

library(rugarch)

library(xts)

library(ADGofTest)

library(qqtest)

library(copula)

library(qrmdata)

library(qrmtools)

library(dplyr)

datax<-read.csv(file="1.csv",sep=",",header=TRUE)

x<-datax$SZI

y<-datax$SH

o<-datax$HSI

plot(x,y)

plot(x,o)

plot(y,o)

data<-data.frame(x,y,o)

S<-data

### 2 Fitting marginal ARMA(1,1)-GARCH(1,1) models with

         standardized t residuals

## Build negative log-returns

X <- -returns(S) # -log-returns

Y <- returns(S)

## Basic plot

plot.zoo(X, main = "Negative log-returns")

plot.zoo(Y, main = "log-returns")

## Fit marginal time series

uspec <- rep(list(ugarchspec(distribution.model = "std")), ncol(X))

fit.ARMA.GARCH <- fit_ARMA_GARCH(X, ugarchspec.list = uspec)

stopifnot(sapply(fit.ARMA.GARCH$error, is.null)) # NULL = no error

if(FALSE)

  fit.ARMA.GARCH$warning

## => Warning comes from finding initial values and can

         be ignored here

fits <- fit.ARMA.GARCH$fit # fitted models

Z <- as.matrix(do.call(merge, lapply(fits, residuals, standardize

        = TRUE))) # grab out standardized residuals

colnames(Z) <- colnames(S)

(nu.mar <- vapply(fits, function(x) x@fit$coef[["shape"]],

        NA_real_)) # vector of estimated df

n <- nrow(X) # sample size

d <- ncol(X) # dimension

### Fitting copulas ##########################################

## Compute pseudo-observations from the standardized t residuals

U <- pobs(Z)

pairs2(U, cex = 0.4, col = adjustcolor("black", alpha.f = 0.5))

## Compare various bivariate copulas

fit.nc <- fitCopula(normalCopula(dim= d),  data = U)

fit.tc <- fitCopula(tCopula(dim= d),       data = U) # df of

        freedom are estimated, too

fit.cl <- fitCopula(claytonCopula(dim= d), data = U)

fit.gc <- fitCopula(gumbelCopula(dim=d),  data = U)

## Comparing the likelihoods

sort(c(nc = fit.nc@loglik, tc = fit.tc@loglik, cl = fit.cl@loglik

        , gc = fit.gc@loglik),

     decreasing = TRUE)

### 3 Fitting copulas #########################################

## Fitting a Gumbel copula

fit.gc <- fitCopula(gumbelCopula(dim = d),

                    data = U, method = "mpl")

fit.gc@estimate # estimated copula parameter

gc <- fit.gc@copula # fitted copula

## Compute matrices of pairwise Kendall's tau and upper

        tail-dependence coefficients

p2P(tau(gc), d = d)

p2P(lambda(gc)["upper"], d = d)

## Fitting a Clayton copula

fit.cl <- fitCopula(claytonCopula(dim = d),

                    data = U)

fit.cl@estimate # estimated copula parameter

cl <- fit.cl@copula # fitted copula

## Compute matrices of pairwise Kendall's tau and lower

        tail-dependence coefficients

p2P(tau(cl), d = d)

p2P(lambda(cl)["lower"], d = d)

## Fitting a Normal copula

fit.nc <- fitCopula(normalCopula(dim = d),

                    data = U)

fit.nc@estimate # estimated copula parameter

nc <- fit.nc@copula # fitted copula

## Compute matrices of pairwise Kendall's tau and lower

        tail-dependence coefficients

p2P(tau(nc), d = d)

p2P(lambda(nc), d = d)

## Fitting a t copula

fit.tc <- fitCopula(tCopula(dim = d, dispstr = "un"),

                    data = U, method = "itau.mpl")

(nu <- tail(fit.tc@estimate, n = 1))

# estimated degrees of freedom nu

(P <- p2P(head(fit.tc@estimate, n = -1)))

# estimated correlation matrix

tc <- fit.tc@copula # fitted copula

## Compute matrices of pairwise Kendall's tau and upper

        tail-dependence coefficients

p2P(tau(tc))

p2P(lambda(tc)[(choose(d,2)+1):(d*(d-1))])

### 4 Goodness-of-fit ##########################################

## We use the parametric bootstrap here, based on the

        maximum pseudo-likelihood

## estimator.

N <- 200

## Check the Gumbel copula

gof.gc <- gofCopula(gc, x = U, N = N) # warning also because

        the copula doesn't fit well here

stopifnot(gof.gc$p.value < 0.05) # if occur=> rejection

## Check the Clayton copula

gof.cl <- gofCopula(cl, x = U, N = N) # warning also because

        the copula doesn't fit well here

stopifnot(gof.cl$p.value < 0.05) # if occur=> rejection

## Check the normal copula

gof.nc <- gofCopula(nc, x = U, N = N) # warning also because

        the copula doesn't fit well here

stopifnot(gof.nc$p.value < 0.05) # if occur=> rejection

## Check the t copula is different from the three above since

        t copula contains two integers.

U.Rsnbl_tc <- cCopula(U, copula = tc)

pairs2(U.Rsnbl_tc, cex = 0.4, col = adjustcolor("black",

        alpha.f = 0.5),main ="t copula") # looks ok

## normal copula also checked

U.Rsnbl_nc <- cCopula(U, copula = nc)

pairs2(U.Rsnbl_nc, cex = 0.4, col = adjustcolor("black",

        alpha.f = 0.5),main ="normal copula") # looks not good

## clayton copula also checked

U.Rsnbl_cl <- cCopula(U, copula = cl)

pairs2(U.Rsnbl_cl, cex = 0.4, col = adjustcolor("black",

        alpha.f = 0.5),main="clayton copula") # looks not good

## Gumbel copula also checked

U.Rsnbl_gc <- cCopula(U, copula = gc)

pairs2(U.Rsnbl_gc, cex = 0.4, col = adjustcolor("black",

        alpha.f = 0.5),main ="gumbel copula") # looks ok

## Map it to a K_d distribution ("kay") and do a Q-Q plot

U.Rsnbl.K <- sqrt(rowMeans(qnorm(U.Rsnbl_tc)^2))

# map to a kay distribution

pK <- function(q, d) pchisq(d*q*q, df = d)

# df of a K_d distribution

AD <- ad.test(U.Rsnbl.K, distr.fun = pK, d = d)

# compute an AD test for t copula

stopifnot(AD$p.value >= 0.05)# contrary to the above test

## A (sophisticated) Q-Q plot, Attetion: we have to clear the

        chart space since we defined above.

qqtest(U.Rsnbl.K, dist = "kay", df = d, nreps = 1000, pch = 1,

       col = adjustcolor("black", alpha.f = 0.5), main = "good",

       xlab = substitute(K[dof]~"quantiles", list(dof = d)))

# => looks ok

\end{lstlisting}

\subsection{t\ copula的蒙特卡罗}

\begin{lstlisting}

## t copula montecarlo

library(mvtnorm)

tcsample<-rmvt(2000,sigma=Sigma,df=v,type=c("shifted","Kshirsagar")[2])

tcsample[,2]<-

(tcsample[,2]-min(tcsample[,2]))/(max(tcsample[,2])-min(tcsample[,2]))

tcsample[,3]<-

(tcsample[,3]-min(tcsample[,3]))/(max(tcsample[,3])-min(tcsample[,3]))

tcsample[,1]<-

(tcsample[,1]-min(tcsample[,1]))/(max(tcsample[,1])-min(tcsample[,1]))

tcsample1<-(apply(tcsample,1,sum))

quantile(tcsample1, probs = seq(0,1, 0.01), type = 7)

\end{lstlisting}

% vim:ts=4:sw=4

举报

相关推荐

0 条评论