temp/iniziale.R

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
lampoverde/Der documentation built on Jan. 8, 2018, 12:01 p.m.