# temp/iniziale.R In lampoverde/Der: Derivatives and Structured Products

```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)){
lapply(funz, function(f) f(pay*disc))
}

perf<-function(x0, xT){xT/x0}
rend<-function(x0, xT){xT/x0-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 )
#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))
}

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
}

#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)
}

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));

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)
cat("Prezzo: ", price\$mean, "Std. Err", price\$sd, "\n")

tmp<-apply(sim, c(1,2), 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
```
lampoverde/Der documentation built on Jan. 8, 2018, 12:01 p.m.