Coursework的代码

install.packages("quantmod")

library(quantmod)

library(xts)

library(zoo)

start<-as.Date("2018-01-01")

end<-as.Date("2019-02-19")

getSymbols("JPM",src="yahoo",from=start,to=end)

JPM

jpm<-JPM$JPM.Adjusted

head(jpm)

jpm.sr<-diff(jpm)/lag(jpm,k=-1)

jpm.sr<-jpm.sr[c(2:282)]

jpm.sr

library(psych)

describe(jpm.sr)

summary(jpm.sr)

sd(jpm.sr,na.rm=T)

str(jpm.sr)

plot(jpm.sr)

#Basic Historical Simulation

Pt<-1000000

loss<--Pt*jpm.sr

head(loss)

v95<-quantile(loss,0.95)

v95

t.loss<-loss[loss>v95]

ES95<-mean(t.loss)

ES95

#Bootstrapping

install.packages("boot")

library(boot)

var95.boot <- function(data,index){

  var.boot <- quantile(data[index],0.95)

  return(var.boot)

}

es95.1 <- function(data,index){

  boot.loss <- data[index]

  boot.var <- quantile (boot.loss, 0.95)

  es.loss <- boot.loss[boot.loss>boot.var]

  boot.es <- mean(es.loss)

  return(boot.es)

}

b.v<-boot(loss,var95.boot,R=1000)

b.v

plot(b.v)

boot.ci(b.v,type="basic")

b.es<-boot(loss,es95.1,R=1000)

b.es

#AW-HS

HS_aw <- function(jpm.sr, P = 1000000, lam = 0.98, ci = 0.95){

  alpha <- 1-ci

  n <- length(jpm.sr)

  df <- data.frame(age = seq(n,1), jpm.sr, loss = -P*jpm.sr)

  df$w <- with(df,((1-lam)*lam^(age-1))/(1-lam^n))

  df <- df[order(df$loss, decreasing = TRUE),]

  df$cumw <- with(df, cumsum(w))

  VaR <- df[which.min(abs(df$cumw - alpha)),'loss']

  df2 <- df[df$loss > VaR,]

  df2$nw <- with(df2, w/sum(w))

  ES <- with(df2, sum(nw * loss))

  res.vec <- c(VaR = VaR, ES = ES)

  names(res.vec) <- paste0(names(res.vec),100*ci)

  return(res.vec)

}

HS_aw(as.vector(jpm.sr))

#Hull-White

library(fGarch)

jpm.dm<-jpm.sr-mean(jpm.sr)

gf<-garchFit(data = jpm.dm)

vol<-volatility(gf)

vol.p<-vol[length(vol)]

st<-vol.p/vol

jpm.adj<-jpm.dm*st

jpm.adj<-jpm.adj+mean(jpm.sr)

jpm.adj

Pt<-1000000

loss1<--Pt*jpm.adj

v95.2<-quantile(loss,0.95)

v95.2

t2.loss<-loss[loss>v95.2]

ES95.2<-mean(t2.loss)

ES95.2

library(boot)

b.v.2<-boot(loss,var95.boot,R=1000)

b.v.2

plot(b.v.2)

boot.ci(b.v.2,type = "basic")

b.es<-boot(loss,es.boot,R=1000)

b.es

plot(b.es)

boot.ci(b.es,type = "basic")

#a normal distribution without voatility adjustment

library(ghyp)

mu<-mean(jpm.sr)

mu

sigma<-sd(jpm.sr)

sigma

g<-gauss(mu,sigma)

g

mean(g)

vcov(g)

plot(g,type="l")

qghyp(0.95,g)

Pt<-1000000

g.l<-transform(g,multiplier=-Pt)

g.l

plot(g.l,type='l')

ESghyp(0.95,g.l,distr="loss")

qghyp(0.95,g.l)

# a normal distribution with voatility adjustment

mu.adj<-mean(jpm.adj)

mu.adj

sigma.adj<-sd(jpm.adj)

sigma.adj

g.adj<-gauss(mu.adj,sigma.adj)

g.adj

mean(g.adj)

vcov(g.adj)

plot(g.adj,type="l")

qghyp(0.95,g.adj)

Pt<-1000000

g.adj.l<-transform(g.adj,multiplier=-Pt)

g.adj.l

plot(g.adj.l,type='l') 

ESghyp(0.95,g.adj.l,distr = "loss")

qghyp(0.95,g.adj.l)

#an appropriate distribution without volatility adjustment

stepAIC.ghyp(jpm.sr)

gf<-fit.hypuv(jpm.sr,symmetric = TRUE)

gf

plot(gf)

hist(gf)

qqghyp(gf,gaussian = FALSE)

s1<-rghyp(10000,gf)

ks.test(s1,jpm.sr)

port1<-1

loss.distribution<-transform(gf,multiplier=-port1)

loss.distribution

ESghyp(0.95,loss.distribution,distr = "loss")

qghyp(0.95,loss.distribution)

# an appropriate distribution with volatility adjustment

stepAIC.ghyp(jpm.adj)

tf<-fit.NIGuv(jpm.adj,symmetric = TRUE)

tf

plot(tf)

hist(tf)

qqghyp(tf,gaussian = FALSE)

s2<-rghyp(10000,tf)

ks.test(s2,jpm.adj)

Pt<-1

loss.distribution.adj<-transform(tf,multiplier=-Pt)

loss.distribution.adj

ESghyp(0.95,loss.distribution.adj,distr = "loss")

qghyp(0.95,loss.distribution.adj)

##two assets

library(quantmod)

library(xts)

library(zoo)

start<-as.Date("2018-01-01")

end<-as.Date("2019-02-19")

getSymbols(c("JPM","AAPL"),src ="yahoo",from=start,to=end)

all<-cbind(JPM[,"JPM.Adjusted"],AAPL[,"AAPL.Adjusted"])

head(all)

##calculate returns

all.ret<-diff(all)/lag(all,k=-1)

all.ret<-all.ret[c(2:282)]

head(all.ret)

summary(all.ret)

plot(all.ret,col=c("blue","red"),legend.loc="top")

library(fBasics)

basicStats(all.ret)

cor(all.ret,use = "complete.obs")

all.ret.df<-as.data.frame(all.ret)

boxplot(all.ret.df)

all.ret

##basic historical simulation

P<-2000000

w<-c(0.5,0.5)

Pw<- -P*w

loss.all<-rowSums(t(Pw*t(all.ret)))

var953<-quantile(loss.all,0.95)

var953

t.loss<-loss.all[loss.all>var953]

ES953<-mean(t.loss)

ES953

##bootstrapping

library(boot)

var.boot <- function(data, index){

  t <- quantile(data[index], 0.95)

  return(t)

}

b.var.all1<- boot(loss.all, var.boot, R=1000)

b.var.all1

plot(b.var.all1)

boot.ci(b.var.all1,type="basic")

es.boot <- function(data,index){

  loss.samp <- data[index]

  var <- quantile(loss.samp,0.95)

  es <- mean(loss.samp[loss.samp>var])

  return(es)

}

b.es.all1 <- boot(loss.all,es.boot,R=1000)

b.es.all1

plot(b.es.all1)

boot.ci(b.es.all1,type="basic")

##Age-weighted historical simulation

HS_aw <- function(r, loss, lam = 0.98, ci = 0.95){

  alpha <- 1-ci

  n <- length(r)

  df <- data.frame(age = seq(n,1), r, loss )

  df$w <- with(df,((1-lam)*lam^(age-1))/(1-lam^n))

  df <- df[order(df$loss, decreasing = TRUE),]

  df$cumw <- with(df, cumsum(w))

  VaR <- df[which.min(abs(df$cumw - alpha)),'loss']

  df2 <- df[df$loss > VaR,]

  df2$nw <- with(df2, w/sum(w))

  ES <- with(df2, sum(nw * loss))

  res.vec <- c(VaR = VaR, ES = ES)

  names(res.vec) <- paste0(names(res.vec),100*ci)

  return(res.vec)

}

HS_aw(all.ret,loss.all)

##Hull-white historical simulation

aapl<-AAPL$AAPL.Adjusted

head(aapl)

aapl.sr<-diff(aapl)/lag(aapl,k=-1)

aapl.sr<-aapl.sr[c(2:282)]

aapl.sr

library(fGarch)

aapl.sr.dm<-aapl.sr-mean(aapl.sr)

gf.aapl<- garchFit(data=aapl.sr.dm)

vol.aapl<-volatility(gf.aapl)

vol.p.aapl<-vol.aapl[length(vol.aapl)]

st.aapl<-vol.p.aapl/vol.aapl

aapl.adj<-aapl.sr.dm*st.aapl

aapl.adj<-aapl.adj+mean(aapl.sr)

aapl.jpm.adj<-cbind(aapl.adj,jpm.adj)

##calculate losses and risk measures

P<-2000000

w<-c(0.5,0.5)

Pw<- -P*w

Loss.adj<-rowSums(t(Pw*t(aapl.jpm.adj)))

Var2.all<-quantile(Loss.adj,0.95)

Var2.all

t.Loss.adj<-Loss.adj[Loss.adj>Var2.all]

ES3.all<-mean(t.Loss.adj)

ES3.all

##bootstrapping

b.var.all1<- boot(Loss.adj,var.boot, R=1000)

b.var.all1

plot(b.var.all1)

boot.ci(b.var.all1,type="basic")

b.es.all1 <- boot(Loss.adj,es.boot,R=1000)

b.es.all1

plot(b.es.all1)

boot.ci(b.es.all1,type="basic")

#Parametric,using the Normal distribution without volatility adjustment

Mean<-colMeans(all.ret)

Mean

Cov.all<- cov(all.ret)

Cov.all

G<-gauss(mu=Mean,sigma=Cov.all)

G

P<-2000000

wa<-0.5

wb<-0.5

mult<- -P*c(wa,wb)

mult

G.l<-transform(G,multiplier=mult)

G.l

ESghyp(0.95,G.l,distr = "loss")

qghyp(0.95,G.l)

##Parametric,using the Normal distribution with volatility adjustment

library(ghyp)

Mean.adj<-colMeans(aapl.jpm.adj)

Mean.adj

Cov.all.adj<-cov(aapl.jpm.adj)

Cov.all.adj

G.adj<-gauss(mu=Mean.adj,sigma=Cov.all.adj)

G.adj

P<-2000000

wa<-0.5

wb<-0.5

mult<- -P*c(wa,wb)

mult

G.l.adj<-transform(G.adj,multiplier=mult)

G.l.adj

ESghyp(0.95,G.l.adj,distr = "loss")

qghyp(0.95,G.l.adj)

##Parametric,using an appropriate distribution without volatility adjustment

stepAIC.ghyp(all.ret)

mf<-fit.NIGmv(all.ret)

mf

#select the variable I want the plot of,

plot(mf[1],type="l")

plot(mf[2],type="l")

hist(mf[1])

qqghyp(mf[1])

hist(mf[2])

qqghyp(mf[2])

#bivariate plot

pairs(mf)

#use Peacock test to check

#check

library(Peacock.test)

S1<-rghyp(10000,mf)

peacock2(all.ret,S1)

#The model is adequate (at 5% significance level) if the test statistic is less than than 1.83, for effective sample size greater than 50.

#calculate effective sample size

n1 <- dim(all.ret)[1]

n2 <- dim(S1)[1]

eff <- (n1*n2)/(n1+n2)

eff

#risk measures

P1 <- 2

wa <- 0.5

wb <- 0.5

mult <- -P1*c(wa,wb)

ld.all<-transform(mf,multiplier=mult)

ld.all

ESghyp(0.95,ld.all,distr = "loss")

qghyp(0.95,ld.all)

##Parametric, using an appropriate distribution with volatility adjustment

stepAIC.ghyp(aapl.jpm.adj)

mf.adj<-fit.tmv(aapl.jpm.adj)

mf.adj

#select the variable I want the plot of,

plot(mf.adj[1],type="l")

plot(mf.adj[2],type="l")

hist(mf.adj[1])

qqghyp(mf.adj[1])

hist(mf.adj[2])

qqghyp(mf.adj[2])

#bivariate plot

pairs(mf.adj)

#use Peacock test to check

library(Peacock.test)

S2<-rghyp(10000,mf.adj)

peacock2(aapl.jpm.adj,S2)

#The model is adequate (at 5% significance level) if the test statistic is less than than 1.83, for effective sample size greater than 50.

#calculate effective sample size

n1.adj <- dim(aapl.jpm.adj)[1]

n2.adj <- dim(S2)[1]

eff.adj <- (n1*n2)/(n1+n2)

eff.adj

#risk measures

P1 <- 2

wa <- 0.5

wb <- 0.5

mult <- -P1*c(wa,wb)

ld.all.adj<-transform(mf.adj,multiplier=mult)

ld.all.adj

ESghyp(0.95,ld.all.adj,distr = "loss")

qghyp(0.95,ld.all.adj)

###(c) find ES in the two-asset case,using basic HS, allowing the weights to vary

#using Hull-White

hsj <- function(wa, P, rets){

  w <- c(wa, 1-wa)

  Pw <- -P*w

  loss <- rowSums(t(Pw * t(rets)))

  var <- quantile(loss, 0.95)

  es <- mean(loss[loss>var])

  return(es)

}

#we want ES for weights from 0-1 at intervals of 0.01

weights<- seq(0, 1, by = 0.01)

sapply(weights, hsj,2000000,aapl.jpm.adj)

ES.weights <- sapply(weights, hsj, 2000000, aapl.jpm.adj)

ES.weights

plot(weights,ES.weights, type = 'l')

#find the minimum ES(risk exposure)

which.min(ES.weights)

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 213,417评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,921评论 3赞 387
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 158,850评论 0赞 349
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,945评论 1赞 285
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,069评论 6赞 385
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,188评论 1赞 291
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,239评论 3赞 412
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,994评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,409评论 1赞 304
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,735评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,898评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,578评论 4赞 336
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,205评论 3赞 317
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,916评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,156评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,722评论 2赞 363
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,781评论 2赞 351

推荐阅读更多精彩内容