library(mvtnorm)
library(plyr)
gbm<-function(T, r, sd, s0=1, step=1){
# Genera un path (GBM vol costante)
# Far restituire anche t
v<-cumprod(c(s0, exp((r - sd^2/2) / (252/step) + sd*rnorm(round(T*252/step,0)) / sqrt (252/step))))
t<-seq(0,T,by=step/252)
list(v=v,t=t)
}
gbmVector<-function(T, r, sd, correl, s0=rep(1,nrow(sigma)), step=1){
# Genera un path (GBM vol costante)
# Far restituire anche t
apply(
rbind(s0, exp((r - sd^2/2) / (252/step) + rmvnorm(round(T*252/step,0), sigma=diag(sd)%*%correl%*%diag(sd)) / sqrt (252/step))),
2, cumprod)
}
# FUNZIONI UTILI NEI PAYOFF
pvPayoffT<-function(pay, disc, funz=list(mean=mean, sd=sd)){
# Sconta payoff a scadenza
lapply(funz, function(f) f(pay*disc))
}
perf<-function(x0, xT){xT/x0}
rend<-function(x0, xT){xT/x0-1}
logret<-function(x){log(tail(x,-1)/head(x,-1))}
basket<-function(x,dim=c(1,3), f, ...){apply(x, dim, f, ...)}
# Fanno perdere 10% tempo (usare direttamente head(x,1) è come usare first e last )
#first<-function(x) head(x, 1)
#last<-function(x) tail(x, 1)
# Questo è come usare direttamente x[1]...
first<-function(x) x[1]
last<-function(x) x[length(x)]
# TODO
replicateApply<-function(){
}
# VARI PAYOFF
# Path independent
call<-function(s,strike){max(s-strike, 0)}
put<-function(s,strike){max(strike-s, 0)}
digitalCall<-function(s, strike){ifelse(s>=strike,1,0)}
digitalPut<-function(s, strike){ifelse(s<=strike,1,0)}
barrOption<-function(s, cond, payoff, reb=0){
# cond: it TRUE option active, else rebate
# Example (geared PUT amer and eur barrier)
# barrOption(s=sim[,1], cond=s<0.9, 2*put(sim[nrow(sim),1],1), reb=0)
# barrOption(s=sim[nrow(sim),1], cond=s<0.95, 2*put(sim[nrow(sim),1],1), reb=0)
ifelse(any(eval(substitute(cond))), payoff, reb)
}
capGar<-function(s, gar, gear){
#gar + s[1]*gear*max(0, rend(s[1], s[length(s)]))
gar + s[1]*gear*max(0, rend(first(s), last(s)))
}
disc<-function(t,r){exp(-r*t)}
eval_mc<-function(N, fgen, pargen, fpay, parpay,fpath=function(x){x},fdisc=function(x){1},pardisc){
res<-val<-array(NA, N)
for(i in 1:N){
path<-do.call(fgen, pargen)
parpay[[1]]<-fpath(path$v)
res[i]<-do.call(fpay, parpay)
pardisc[[1]]<-fpath(path$t)
val[i]<-do.call(fdisc,pardisc)*res[i]
}
list(m=mean(val),s=sd(val))
}
callSPread<-function(s,strike1,strike2){call(s,strike1)-call(s,strike2)}
putSPread<-function(s,strike1,strike2){put(s,strike1)-put(s,strike2)}
tmp2<-Vectorize(FUN = callSPread)
tmp2<-function(x){Vectorize(FUN = callSPread)(s=x,1,1.2)}
tmp<-function(x)tmp2(s=x,1)
windows()
curve(tmp2,from=.5,to=1.7,lwd=2,col='red'); grid()
y<-gbm(T=1, r=0.01, sd=0.2, s0=1, step=1/252)
# =========================================================================================
disc<-function(t,r){exp(-r*t)}
discFix<-function(r){function(t) exp(-r*t)}
gbmTimeVar<-function(T, f_disc=function(t) 1, f_dvd=function(t) 0, f_sd=function(t) 0,
f_fxvol=function(t) 0, fxcor=0, s0=1, step=1,f_rand=rnorm,p_rand=list(n=NULL)){
# Genera un path (GBM r e vol time varying)
# rendere custom il path (punti aggiuntivi etc)
t<-seq(0,T,by=step/252)
N<-length(t)
dt<-t[-1]-t[-N]
sd<-sapply(t[-1],f_sd)
dvd<-sapply(t[-1],f_dvd)
disc<-sapply(t, f_disc)
fxvol<-sapply(t[-1], f_fxvol)
p_rand[[1]]<-N-1
if(is.null(f_rand)){
function(z){
tmp<-cumprod(c(s0, disc[-N]/disc[-1]*exp( (-sd^2/2 - dvd - fxcor*fxvol*sd)*dt + sd*z*sqrt(dt))))
}
}
else{
function(){
tmp<-cumprod(c(s0, disc[-N]/disc[-1]*exp( (-sd^2/2 - dvd - fxcor*fxvol*sd)*dt + sd*do.call(f_rand,p_rand)*sqrt(dt))))
}
}
}
gbmVectTimeVar<-function(T, gens, step=1,f_rand=rmvnorm,p_rand=list(n=NA,sigma=diag(length(gens)))){
# rendere custom il path (punti aggiuntivi etc)
t<-seq(0,T,by=step/252)
N<-length(t)
dt<-t[-1]-t[-N]
K<-length(gens)
s0<-array(NA,K)
DISC<-matrix(NA,N,K)
DVD<-FXVOL<-SD<-matrix(NA,N-1,K)
i<-1
for(g in gens){
s0<-environment(g)$s0
DISC[,i]<-environment(g)$disc
DVD[,i]<-environment(g)$dvd
FXVOL[,i]<-environment(g)$fxvol
SD[,i]<-environment(g)$sd
i<-i+1
}
p_rand[[1]]<-N-1
if(is.null(f_rand)){
function(z){
apply(
rbind(s0, DISC[-N,]/DISC[-1,]*exp( (-SD^2/2 - DVD - fxcor*FXVOL*SD)*dt + SD*z*sqrt(dt))),
2, cumprod)
}
}
else{
function(){
apply(
rbind(s0, DISC[-N,]/DISC[-1,]*exp( (-SD^2/2 - DVD - fxcor*FXVOL*SD)*dt + SD*do.call(f_rand, p_rand)*sqrt(dt))),
2, cumprod)
}
}
}
eval_mc<-function(N, fgen, fpay, parpay,fpath=function(x){x},fdisc=function(x){1}){
t<-environment(fgen)$t
res<-val<-array(NA, N)
for(i in 1:N){
path<-fgen()
parpay[[1]]<-fpath(path)
res[i]<-do.call(fpay, parpay)
tval<-fpath(t)
val[i]<-fdisc(tval)*res[i]
}
list(m=mean(val),s=sd(val))
}
# TODO: con f_rand NULL non funziona
gen1<-gbmTimeVar(1,f_disc=disc,f_dvd=function(t) .0363, f_sd=function(t) .1963)
gen2<-gbmTimeVar(1,f_disc=disc,f_dvd=function(t) .0363, f_sd=function(t) .1963)
gen3<-gbmTimeVar(1,f_disc=disc,f_dvd=function(t) .0363, f_sd=function(t) .1963)
gens<-c(gen1,gen2,gen3)
gensVect<-gbmVectTimeVar(1, gens, step=1,f_rand=rmvnorm,p_rand=list(n=NA,sigma=diag(length(gens))))
Y<-gensVect()
plot.ts(Y, plot.type = 'multiple')
eval_mc_correlated<-function(N, fgens, correl, fpay, parpay,fpath=function(x){x},fdisc=function(x){1}){
t<-environment(fgen)$t
res<-val<-array(NA, N)
for(i in 1:N){
path<-fgen()
parpay[[1]]<-fpath(path)
res[i]<-do.call(fpay, parpay)
tval<-fpath(t)
val[i]<-fdisc(tval)*res[i]
}
list(m=mean(val),s=sd(val))
}
disc<-discFix(-.02)
gen<-gbmTimeVar(1,f_disc=disc,f_dvd=function(t) .0363, f_sd=function(t) .1963)
eval_mc(N=100, gen, call, list(s=NA,strike=1.0),fpath=last, fdisc=disc)
MCiter<-100
y<-matrix(NA,environment(gen)$N,MCiter)
for(i in 1:MCiter){
y[,i]<-gen()
}
plot.ts(y,plot.type = 'single', col='lightblue'); grid()
# =========================================================================================
# Path dependent con payoff a scadenza
certPlus<-function(s, gar,barr){
#ifelse(min(s)>barr,max(gar,s[length(s)]),s[length(s)])
ifelse(min(s)>barr,max(gar,last(s)),last(s)) # Perde 10% tempo
}
# Path dependent su basket, con payoff a scadenza
#certPlus<-function(s, gar,barr,fbarr,fperf, ...){
# ifelse(min(fbarr(s,...))>barr,max(gar,fperf(s[length(s)],...)),fperf(s[length(s)]))
#}
# Path dependent con payoff eventualmente anticipato
# Gestire caso barriera europea
autocallable<-function(sim, t, obs, trig, barr, cpn, cap){
w<-match(obs, t)
wmin<-which.min(sim[w]>=trig)
tmin<-t[w[wmin]]
coup<-wmin*cpn
prot<-ifelse(min(sim)<=barr,sim[length(sim)],cap)
ifelse(coup>0,{out<-c(cap+cpn,tmin,wmin)},{out<-c(prot,t[length(t)],length(obs))})
names(out)<-c('pay','t','i')
return(out)
}
# BASKET ===================================================================
tmp<-aaply(sim, c(1,2), mean)
out<-aaply(sim1, 2, each(mean, sd))
sim1<-replicate(N, gbm(T=3, r=0.03, sd=0.2, s0=1, step=3));
# Basket option
sim<-replicate(N, gbmVector(T=3, r=0.03, sd=c(0.2,.1),correl=matrix(c(1,.7,.7,1),2,2), s0=c(1,1), step=3));
f1<- function(x) call(x, strike=1)
price<-pvPayoffT(apply(basket(sim,2,mean), 3, f), 0.92);
cat("Prezzo: ", price$mean, "Std. Err", price$sd, "\n")
tmp<-apply(sim, c(1,2), mean)
tmp<-basket(x=sim, f=mean)
# ESEMPI ===================================================================
N<-20000
sim<-replicate(N, gbm(T=3, r=0.03, sd=0.2, s0=1, step=3));
pricePlus<-pvPayoffT(apply(sim, 2, certPlus, gar=1,barr=0.92), 0.92);
cat("Prezzo Plus: ", pricePlus$mean, "Std. Err", pricePlus$sd, "\n")
# Cert. plus (20'000 iter, 3 anni, sampling ogni 1gg) ...vediamo quanto ci mette
system.time({
sim<-replicate(N, gbm(T=3, r=0.03, sd=0.2, s0=1, step=3));
pricePlus<-pvPayoffT(apply(sim, 2, certPlus, gar=1,barr=0.92), 0.92);
cat("Prezzo Plus: ", pricePlus$mean, "Std. Err", pricePlus$sd, "\n")
}
)
# Vari payoff (20'000 iter, 3 anni, sampling ogni 1gg) ...scrittura piu' prolissa
system.time({
sim<-replicate(N, gbm(T=3, r=0.03, sd=0.2, s0=1, step=3));
#f<- function(x) certPlus(x, gar=1,barr=0.92)
f1<- function(x) call(x, strike=1)
f2<- function(x) call(x, strike=1)
f <- function(x) f1(x) - f2(x)
f <- function(x) list(f1(x), f2(x), f1(x)+f2(x))
#f<- function(x) barrOption(s=x, cond=s<0.9, 2*put(x[length(x)],strike=1), reb=0)
price<-pvPayoffT(apply(sim, 2, f), 0.92);
cat("Prezzo: ", price$mean, "Std. Err", price$sd, "\n")
}
)
# Grafico 100 scenari
plot.ts(sim[,sample(1:N,100)], plot.type='single', col='lightblue')
abline(a=0.9,b=0, lty=2,col='red')
grid()
# Usiamo gli stessi scenari per un capitale garantito
payoffCG<-apply(sim, 2, capGar,gar=1,gear=0.25)
priceCG<-pvPayoffT(payoffCG, 0.92);
cat("Prezzo CG: ", priceCG$mean, "\n")
# Autocallable (df e obs includono la scadenza)
dfpay<-c(0.985,0.965,0.945,0.925, 0.905,0.885)
simt<-seq(0,3,length.out=nrow(sim)) #Far restituire dal simulatore...
payoffAC<-apply(sim, 2, autocallable, simt, obs=c(0.5,1,1.5,2,2.5,3), trig=1, barr=0.8, cpn=0.05, cap=1);
priceAC<-pvPayoffT(payoffAC['pay',], dfpay[payoffAC['i',]]);
cat("Prezzo AC: ", priceAC$mean, "\n")
pvPayoffT(sapply(sim[nrow(sim),],call,strike=1.05), 0.92);
pvPayoffT(sapply(sim[nrow(sim),],put,strike=0.95), 0.92);
#FINIRE
payoff.plot<-function(f,s,m,M,N=100, ...){
ftmp<-Vectorize(f,s)
s<-seq(m,M,length.out=N)
v<-ftmp(s)
plot(s,v, type='l', lwd=2)
grid()
}
#
function S = mod_pathgener(tempo,K,Nasset,x,Scurr,r,d,sigma,sigma_fx,rho_fx)
% help
%
% in caso di local vol, "sigma" è la first vol
dd=365;
% S=zeros(Nasset,K);
% drift=exp((r-d)*tempo(1)/dd);
% S(:,1)=Scurr.*drift.*exp(((-rho_fx.*sigma.*sigma_fx)*(tempo(1)/dd)-0.5*(tempo(1)/dd)*sigma.^2)+sqrt(tempo(1)/dd)*sigma.*x(:,1));
% for j=2:K
% dt=tempo(j)-tempo(j-1);
% drift=exp((r-d)*dt/dd);
% S(:,j)=S(:,j-1).*drift.*exp(((-rho_fx.*sigma.*sigma_fx)*(dt/dd)-0.5*(dt/dd)*sigma.^2)+sqrt(dt/dd)*sigma.*x(:,j));
% end
moto=zeros(Nasset,K);
moto(:,1)=Scurr.*exp((r-d)*tempo(1)/dd).*exp(((-rho_fx.*sigma.*sigma_fx)*(tempo(1)/dd)-0.5*(tempo(1)/dd)*sigma.^2)+sqrt(tempo(1)/dd)*sigma.*x(:,1));
dt=tempo(2:end)-tempo(1:end-1);
moto(:,2:end)=exp((r-d-rho_fx.*sigma.*sigma_fx-0.5*sigma.^2)*dt/dd+diag(sigma)*x(:,2:end).*sqrt(repmat(dt,Nasset,1)/dd));
S=cumprod(moto,2);
function S = mod_pathgener_fc(tempo,K,Nasset,x,Scurr,DFx,DFy,sigma,sigma_fx,rho_fx)
% help
%
% in caso di local vol, "sigma" è la first vol
dd=365;
DFy_interp=interp1(DFx,DFy,tempo/dd,'linear');
drift=DFy_interp(1:end-1)./DFy_interp(2:end);
S=zeros(Nasset,K);
S(:,1)=Scurr.*(1/DFy_interp(1)).*exp(((-rho_fx.*sigma.*sigma_fx)*(tempo(1)/dd)-0.5*(tempo(1)/dd)*sigma.^2)+sqrt(tempo(1)/dd)*sigma.*x(:,1));
for j=2:K
dt=tempo(j)-tempo(j-1);
S(:,j)=S(:,j-1).*drift(j-1).*exp(((-rho_fx.*sigma.*sigma_fx)*(dt/dd)-0.5*(dt/dd)*sigma.^2)+sqrt(dt/dd)*sigma.*x(:,j));
end
function [S, lv] = mod_pathgener_lv(tempo,K,Nasset,x,Scurr,r,d,sigma,X_fine_all,Y_fine_all,local_fine_all,sigma_fx,rho_fx,dove)
% help
%
% in caso di local vol, "sigma" è la first vol
dd=365;
lv=zeros(Nasset,K);
lv(:,1) = sigma;
S=zeros(Nasset,K);
drift=exp((r-d)*tempo(1)/dd);
S(:,1)=Scurr.*drift.*exp(((-rho_fx.*sigma.*sigma_fx)*(tempo(1)/dd)-0.5*(tempo(1)/dd)*sigma.^2)+sqrt(tempo(1)/dd)*sigma.*x(:,1));
for j=2:K
dt=tempo(j)-tempo(j-1);
drift=exp((r-d)*dt/dd);
for k=dove
sigma(k)=vicino(X_fine_all{k},Y_fine_all{k},local_fine_all{k},tempo(j-1)/dd,S(k,j-1)/Scurr(k));
end
S(:,j)=S(:,j-1).*drift.*exp(((-rho_fx.*sigma.*sigma_fx)*(dt/dd)-0.5*(dt/dd)*sigma.^2)+sqrt(dt/dd)*sigma.*x(:,j));
lv(:,j) = sigma;
end
function S = mod_pathgener_lv_fc(tempo,K,Nasset,x,Scurr,DFx,DFy,sigma,X_fine_all,Y_fine_all,local_fine_all,sigma_fx,rho_fx,dove)
% help
%
% per il momento, la funzione è valida solo nel caso di sottostante univariato
% "sigma" è la first vol
dd=365;
DFy_interp=interp1(DFx,DFy,tempo/dd,'linear');
drift=DFy_interp(1:end-1)./DFy_interp(2:end);
S=zeros(Nasset,K);
S(:,1)=Scurr.*(1/DFy_interp(1)).*exp(((-rho_fx.*sigma.*sigma_fx)*(tempo(1)/dd)-0.5*(tempo(1)/dd)*sigma.^2)+sqrt(tempo(1)/dd)*sigma.*x(:,1));
for j=2:K
dt=tempo(j)-tempo(j-1);
for k=dove
sigma(k)=vicino(X_fine_all{k},Y_fine_all{k},local_fine_all{k},tempo(j-1)/dd,S(k,j-1)/Scurr(k));
end
S(:,j)=S(:,j-1).*drift(j-1).*exp(((-rho_fx.*sigma.*sigma_fx)*(dt/dd)-0.5*(dt/dd)*sigma.^2)+sqrt(dt/dd)*sigma.*x(:,j));
end
function S = mod_pathgener_lv_fc_multiv(tempo,K,Nasset,x,Scurr,DFx,DFy,sigma,X_fine_all,Y_fine_all,local_fine_all,sigma_fx,rho_fx,dove)
% help
%
% per il momento, la funzione è valida solo nel caso di sottostante univariato
% "sigma" è la first vol
dd=365;
DFy_interp=interp1(DFx,DFy,tempo/dd,'linear');
drift=DFy_interp(1:end-1)./DFy_interp(2:end);
S=zeros(Nasset,K);
S(:,1)=Scurr.*(1/DFy_interp(1)).*exp(((-rho_fx.*sigma.*sigma_fx)*(tempo(1)/dd)-0.5*(tempo(1)/dd)*sigma.^2)+sqrt(tempo(1)/dd)*sigma.*x(:,1));
for j=2:K
dt=tempo(j)-tempo(j-1);
for k=dove
sigma(k)=vicino(X_fine_all{k},Y_fine_all{k},local_fine_all{k},tempo(j-1)/dd,S(k,j-1)/Scurr(k));
end
S(:,j)=S(:,j-1).*drift(j-1).*exp(((-rho_fx.*sigma.*sigma_fx)*(dt/dd)-0.5*(dt/dd)*sigma.^2)+sqrt(dt/dd)*sigma.*x(:,j));
end
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.