-
Notifications
You must be signed in to change notification settings - Fork 0
/
logistic2009_vol_changed.R
214 lines (157 loc) · 6.02 KB
/
logistic2009_vol_changed.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
# Numerical example 1
# This code produce the result in Table 6: The bond price movement against volatility terms
# You need bond coupons from hull white model
library(COVID19)
library(demodelr)
library(matrixStats)
library(ggplot2)
library(forecast)
library(tseries)
library(RQuantLib)
library(quantmod)
library(dplyr)
library(tidyverse)
library(zoo)
library(Sim.DiffProc)
library(rstudioapi)
library(xtable)
graphics.off() # clear all graphs
rm(list = ls()) # remove all files from your workspace
# change directory to where the script located
my_d <- dirname(rstudioapi::getActiveDocumentContext()$path)
setwd(my_d)
#######################################changing volatility########################################
set.seed(123)
Total_infected<-833*10^6
Total_death<-0.4*10^6
g_I<-0.04
sigma<-seq(0.02,0.4,0.02)
Bond_Price<-c()
P_C<-c()
for(j in 1:length(sigma)){
print(paste0("Simulating Scenario ", j , " out of ",length(sigma)))
y1_max<- Total_infected #maximum infected
y2_max<-Total_death# maximum deaths
y1_min<-1 # initial number of infected
y2_min<-1 #initial number of death
numsim<-5000 # number of simulations
n<-1104 # number of days between July 07, 2017 -July 15, 2020 (including start date)
theta_I <-5000
theta_D<-2500
s0_I<-sigma[j] #volatility for infected
r0_I<-g_I #growth rate for infected
K0_I<-y1_max #carrying capacity
x0_I<-y1_min #N0_I
s0_D<-sigma[j] #volatility for death
r0_D<-g_I #growth rate for death
K0_D<-y2_max
x0_D<-y2_min #N0_D
dt<-1 #time step
m<-n #number of steps
M<-numsim #number of paths(simulations)
rho<-0.5 #correlation coefficient between death and infected
#######################################################################
fun_st_logistic<- function(dt,m,rho,M){
fx<-expression(r0_I*x*(1-x/(K0_I )),r0_D*y*(1-y/(K0_D)))
gx<-expression(s0_I*x*(1-x/(K0_I )),s0_D*y*(1-y/(K0_D)))
x0<-c(x0_I,x0_D)
t0<-0
Sigma<-matrix(c(1,rho,rho,1),nrow=2)
set.seed(123)
res<-snssde2d(N=m,M=M,x0=x0,t0=t0,Dt=dt,drift=fx,diffusion=gx,corr=Sigma)
return(res )
}
res<-fun_st_logistic(dt,m,rho,M)
simulated_I<-matrix(unlist(res$X),ncol=numsim,byrow=F)
simulated_D<-matrix(unlist(res$Y),ncol=numsim,byrow=F)
###########################We calculate seven day average for the data set ################
###First find the real increment from previous day to today ##############
Increment_I<-colDiffs(simulated_I, lag=1)
Increment_I<-pmax(Increment_I,0)
nn<-nrow(Increment_I)
index<-1:nn
sim_avg7_I<- matrix(ncol = numsim, nrow = length(7:nn))
sim_GR_I<-matrix(ncol = numsim, nrow = length(7:nn))
for(i in 1:numsim){
obj1<-zoo(Increment_I[,i],index)
seven_avg_I<-rollmean(obj1,7,align='right',fill=0)
seven_sd_I<-rollapply(data = obj1,width=7,FUN=sd,align='right',fill=0)
y11<-as.data.frame(seven_avg_I[7:nn])
y11<-y11[["seven_avg_I[7:nn]"]]
sim_avg7_I[,i]<-y11
y33<-as.data.frame(seven_sd_I[7:nn])
y33<-y33[["seven_sd_I[7:nn]"]]
sim_GR_I[,i]<-(y11-(1.533*y33))
}
Increment_D<-colDiffs(simulated_D, lag=1)
Increment_D<-pmax(Increment_D,0)
sim_avg7_D<- matrix(ncol = numsim, nrow = length(7:nn))
for(i in 1:numsim){
obj1<-zoo(Increment_D[,i],index)
seven_avg_D<-rollmean(obj1,7,align='right',fill=0)
y22<-as.data.frame(seven_avg_D[7:nn])
y22<-y22[["seven_avg_D[7:nn]"]]
sim_avg7_D[,i]<-y22
}
##################################################################
trigger<-0
set_I<-0
set_D<-0
set_GR<-0
common<-0
common_index<-c()
common_time<-c()
for(i in 1:numsim){
set_I<-which(sim_avg7_I[,i]>theta_I)
set_D<-which(sim_avg7_D[,i]>theta_D)
set_GR<-which(sim_GR_I[,i]>0)
common<-Reduce(intersect, list(set_I,set_D,set_GR))
if(length(common)>0){
trigger<-trigger+1
common_index<-append(common_index,i) #in which simulation this happen
common_time<-append(common_time,common[1]) #at what time this happen, then add first time all criteria met
}
common<-NULL
}
#read predicted coupons. last elements is coupon+redemption amount
coupon_data<-readRDS(file = "coupon.rds")
#A function to calculate present value
pv <- function(y, cf, t){ sum(cf/(1+y)^t)}
#y1 is the yield rate calculated early assuming no pandemic.
y1<-0.08673402
BondValue_partial<-c()
for(i in 1:length(common_index)){
trigger_time<-common_time[i]+6 #we add 6, since these are obtained from 7 day average, which set 7->1. So add 6 to bring back to original
coupondf<-subset(coupon_data, coupon_data$Time<=trigger_time)
cf<-as.numeric(coupondf$CF)
t<-coupondf$Time/360
presentV<-pv(y1,cf,t)
BondValue_partial<-append(BondValue_partial,presentV)
}
#E[BondPrice|Trigger]
price1<-mean(BondValue_partial)
#E[BondPrice|No Trigger]
price2<-225 #this is the price advertised by World Bank/paid by investors
#P(H) probability of pandemic
prob1<-0.1393
#P(C|H) probability that trigger activated given pandemic
print(trigger)
prob2<-(trigger/numsim)
# P(C) probability of trigger activated
prob3<-prob1*prob2
#P(C|H') probability that trigger activated when no pandemic is zero.
#P(C') probability that trigger not activated
prob4<-1-prob3
# The price of the bond
#E[Bond Price|Trigger]P (Trigger)+E[Bond Price| No Trigger]P (No Trigger)
price3<-(price1*prob3)+(price2*prob4)
Bond_Price<-append(Bond_Price,price3)
P_C<-append(P_C,prob3)
}
aa1<-round(Bond_Price,4)
aa2<-round(P_C,4)
df_bond<-as.data.frame(rbind(sigma,aa1,aa2))
rownames(df_bond)<-c("sigma","Bond.Price","Prob.Trigger")
xtable(df_bond[,c("V2","V4","V6","V8","V10","V12","V14","V16","V18","V20")])
saveRDS(df_bond, file = "bond_price_vol_change.rds")
mean(Bond_Price) # this is the price investor should pay