VAR_calculations

#Exercise, using montecarlo simulation calculate the 10 day 99% var for a portfolio of 20 million equally invested in call and put options with the following parameters

call<-c(30,50,0.05,0.5,0.75,.02,.75,10/252,.75-(10/252))
put<-c(35,20,0.08,0.8,1,.02,.75,10/252,1-(10/252))
parameters<-rbind(call,put)
colnames(parameters)<-c("S_0","K","mu","sigma","T","r","p","dt","t-dt")
print(parameters)

##      S_0  K   mu sigma    T    r    p         dt      t-dt
## call  30 50 0.05   0.5 0.75 0.02 0.75 0.03968254 0.7103175
## put   35 20 0.08   0.8 1.00 0.02 0.75 0.03968254 0.9603175

#generating realisations of 2 normal distributions with correlation .75 for the error term

n<-1000
set.seed(420)
rho<-.75
mu1<-0;s1<-1
mu2<-0;s2<-1
mu<-c(mu1,mu2)
#covariance matrix
sigma<-matrix(c(s1^2,s1*s2*rho,s1*s2*rho,s2^2),nrow=2)
rownames(sigma)<-c(1,2)
colnames(sigma)<-c(1,2)
library(ggplot2)
library(MASS)
bvn<-mvrnorm(n,mu,sigma)

#simulating stock price,d1, d2 in 10 days

#define s_t_ij+delta t = st_2ij for i observation and asset j
st_2ji<-function(s0,mu,sigma,dtx,error){
  return(s0*exp((mu-(0.5*sigma^2))*dtx+(sigma*sqrt(dtx)*error)))
}
#call stock prices
st_2a<-st_2ji(s0=parameters[1,1],
              mu=parameters[1,3],
              dtx=parameters[1,8],
              error=bvn[,1],
              sigma=parameters[1,4])
#put stock prices
st_2b<-st_2ji(s0=parameters[2,1],
              mu=parameters[2,3],
              dtx=parameters[2,8],
              error=bvn[,2],
              sigma=parameters[2,4])


#defining d1 function
function_d1<-function(st,k,sigma,dtx){ 
  return((log(st/k)+
            ((parameters[1,6]+(sigma^2/2))*dtx)) /
           (sqrt(dtx)*sigma))}
#calculating d1 
d1_a<-function_d1(st=st_2a,
                  k=parameters[1,2],
                  dtx=parameters[1,9],
                  sigma=parameters[1,4])

d1_b<-function_d1(st=st_2b,
                  parameters[2,2],
                  dtx=parameters[2,9],
                  sigma=parameters[2,4])
#defining d2 function
fucntion_d2<-function(d1,sigma,dtx){
  return(d1-(sigma*sqrt(dtx)))
  }
#calculating d2
d2_a<-fucntion_d2(d1=d1_a,
                  sigma=parameters[1,4],
                  dtx=parameters[1,9])  
d2_b<-fucntion_d2(d1=d1_b,
                  sigma=parameters[2,4],
                  dtx=parameters[2,9])

#calculating call price

function_call<-function(st,d1,d2,k,dtx){
  return(st*pnorm(d1)-(k/exp(parameters[1,6]*dtx)*pnorm(d2)))}
c_a<-function_call(st_2a,
                   d1=d1_a,
                   d2=d2_a,
                   k=parameters[1,2],
                   dtx = parameters[1,9])

#calculating put price

function_put<-function(st,d1,d2,k,dtx){
  return((k/exp(parameters[1,6]*dtx)*pnorm(-d2))-(st*pnorm(-d1)))}
p_b<-function_put(st=st_2b,
                  d1=d1_b,
                  d2=d2_b,
                  k=parameters[2,2],
                  dtx = parameters[2,9])

#calculating portfolio #calculating purchase price of call and put

d1_call<-function_d1(st=parameters[1,1],
                     k=parameters[1,2],
                     dtx = parameters[1,5],
                     sigma=parameters[1,4])
d1_put<-function_d1(st=parameters[2,1],
                    k=parameters[2,2],
                    dtx=parameters[2,5],
                    sigma=parameters[2,4])

d2_call<-fucntion_d2(d1=d1_call,
                     sigma=parameters[1,4],
                     dtx = parameters[1,5])
d2_put<-fucntion_d2(d1=d1_put,
                    sigma=parameters[2,4],
                    dtx=parameters[2,5])

call_price<-function_call(st=parameters[1,1],
                          d1=d1_call,
                          d2=d2_call,
                          k=parameters[1,2],
                          dtx=parameters[1,5])
put_price<-function_put(st=parameters[2,1],
                        d1=d1_put,
                        d2=d2_put,
                        k=parameters[2,2],
                        dtx = parameters[2,5])

call_price

## [1] 1.027953

put_price

## [1] 2.743948

#quantity of call/puts bought with 10 million allocation each (in 000's)
q_c<-10000/call_price
q_p<-10000/put_price
V_t<-(2*10000*exp(parameters[1,6]*10/252))-(q_c*c_a+q_p*p_b)

#sorting data to estimate VAR

v_t<-as.data.frame(sort(V_t,decreasing=F))
names(v_t)<-"sorted"
v_lowest_20<-as.data.frame(v_t$sorted[1:20])
names(v_lowest_20)<-"sorted"
#10 day 99% VAR
v_lowest_20$sorted[10]

## [1] -14447.32

library(ggplot2)
ggplot(v_t,aes(x=sorted))+geom_histogram(binwidth = 1750)

#ES 98%
ggplot(v_lowest_20,aes(x=sorted))+geom_histogram(binwidth = 5000)

sum(v_lowest_20)/length(v_lowest_20)

## [1] -301924.9