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)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末靖避，一起剝皮案震驚了整個(gè)濱河市跋理，隨后出現(xiàn)的幾起案子蒸殿，更是在濱河造成了極大的恐慌筷厘，老刑警劉巖，帶你破解...
    沈念sama閱讀 212,454評論 6贊 493
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件宏所，死亡現(xiàn)場離奇詭異酥艳，居然都是意外死亡，警方通過查閱死者的電腦和手機(jī)爬骤，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,553評論 3贊 385
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門充石，熙熙樓的掌柜王于貴愁眉苦臉地迎上來，“玉大人盖腕，你說我怎么就攤上這事赫冬。” “怎么了溃列？”我有些...
    開封第一講書人閱讀 157,921評論 0贊 348
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵劲厌，是天一觀的道長。 經(jīng)常有香客問我听隐，道長，這世上最難降的妖魔是什么雅任？ 我笑而不...
    開封第一講書人閱讀 56,648評論 1贊 284
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮硼婿，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘寇漫。我一直安慰自己刊殉，他們只是感情好州胳，可當(dāng)我...
    茶點(diǎn)故事閱讀 65,770評論 6贊 386
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著栓撞，像睡著了一般。 火紅的嫁衣襯著肌膚如雪瓤湘。 梳的紋絲不亂的頭發(fā)上，一...
    開封第一講書人閱讀 49,950評論 1贊 291
  • 城市分裂傳說
    那天惜索，我揣著相機(jī)與錄音特笋，去河邊找鬼。 笑死猎物，一個(gè)胖子當(dāng)著我的面吹牛，可吹牛的內(nèi)容都是我干的蔫磨。 我是一名探鬼主播，決...
    沈念sama閱讀 39,090評論 3贊 410
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼堤如，長吁一口氣：“原來是場噩夢啊……” “哼！你這毒婦竟也來了搀罢？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 37,817評論 0贊 268
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤榔至，失蹤者是張志新（化名）和其女友劉穎，沒想到半個(gè)月后唧取，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體，經(jīng)...
    沈念sama閱讀 44,275評論 1贊 303
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡邢享，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 36,592評論 2贊 327
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了骇塘。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片掸犬。...
    茶點(diǎn)故事閱讀 38,724評論 1贊 341
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡湾碎，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出介褥，到底是詐尸還是另有隱情，我是刑警寧澤柔滔，帶...
    沈念sama閱讀 34,409評論 4贊 333
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站睛廊，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏超全。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 40,052評論 3贊 316
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一邓馒、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧光酣，春花似錦、人聲如沸救军。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,815評論 0贊 21
  • 一樁弒父案唱遭，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至胆萧，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間跌穗，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 32,043評論 1贊 266
  • 情欲美人皮
    我被黑心中介騙來泰國打工蚌吸， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留，地道東北人羹唠。 一個(gè)月前我還...
    沈念sama閱讀 46,503評論 2贊 361
  • 代替公主和親
    正文 我出身青樓娄昆，卻偏偏與公主長得像，于是被迫代替她去往敵國和親萌焰。 傳聞我的和親對象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 43,627評論 2贊 350

推薦閱讀更多精彩內(nèi)容