bs.binomial<-function(S,K,sigma,r,T,n,type="call") { # S stock price # K strike # sigma annualized volatility in percents # r interest rate (annual) # T time to expiration (fractions of year) q<-1/2+r*sqrt(T/n)/(2*sigma) a<-1+sigma*sqrt(T/n) C<-0 for (j in 0:n){ snj<-a^j*(2-a)^(n-j)*S if (type=="call") Hn<-max(snj-K,0) else Hn<-max(K-snj,0) #C<-C+factorial(n)/(factorial(n-j)*factorial(j))*q^j*(1-q)^(n-j)*Hn C<-C+choose(n,j)*q^j*(1-q)^(n-j)*Hn } C<-(1+r*T/n)^(-n)*C return(C) } bs.delta<-function(S,K,v,r,t,type="call") { # S stock price # K strike # v volatility in percents # r interest rate (annual) # t time to expiration (fractions of year) h<-(log(S/K)+(r+v^2/2)*t)/(v*sqrt(t)) delta<-pnorm(h) if (type=="put") C<-C-S+K*exp(-r*t) return(delta) } bs.implied<-function(P,S,K,r,t,type="call",vlow=0.01,vupp=0.9,tol=0.01) { # P is the price of the option # S stock price # K strike # r interest rate (annual) # t time to expiration (fractions of year) # returns implied volatility in percents vcurlow<-vlow vcurupp<-vupp vcur<-(vcurlow+vcurupp)/2 diff<-vupp-vlow while (diff>=tol){ cur<-bs(S,K,vcur,r,t,type=type) if (cur>=P) vcurupp<-vcur else vcurlow<-vcur vcur<-(vcurlow+vcurupp)/2 diff<-vcurupp-vcurlow } return(vcur) } bs<-function(S,K,v,r,t,type="call") { # S stock price # K strike # v volatility (annualized) # r interest rate (annual) # t time to expiration (fractions of year) if (K==0) C<-S else{ h<-(log(S/K)+(r+v^2/2)*t)/(v*sqrt(t)) C<-S*pnorm(h)-K*exp(-r*t)*pnorm(h-v*sqrt(t)) } if (type=="put") C<-C-S+K*exp(-r*t) return(C) } bs.recursive.amerput<-function(S,K,sigma,r,dt,n) { # S stock price # K strike # sigma annualized volatility in percents # r interest rate (annual) # dt time to expiration (fractions of year) # n is the number of steps u<-1+sigma*sqrt(dt/n) q<-1/2+r*sqrt(dt/n)/(2*sigma) stock.prices<-matrix(0,n+1,1) for (i in 0:n) stock.prices[i+1]<-u^i*(2-u)^(n-i)*S prev<-pmax(K-stock.prices,0) for (k in (n-1):0){ deriv.prices<-matrix(0,k+1,1) stock.prices<-matrix(0,k+1,1) for (i in 0:k) stock.prices[i+1]<-u^i*(2-u)^(k-i)*S for (i in 1:(k+1)) deriv.prices[i]<- max(K-stock.prices[i],(1+r*dt/n)^(-1)*((1-q)*prev[i]+q*prev[i+1])) prev<-deriv.prices } return(prev[1]) } bs.recursive<-function(S,K,sigma,r,T,n,type="amer.put") { # S stock price # K strike # sigma annualized volatility in percents # r interest rate (annual) # T time to expiration (fractions of year) # n is the number of steps u<-1+sigma*sqrt(T/n) q<-1/2+r*sqrt(T/n)/(2*sigma) stock.prices<-matrix(0,n+1,1) for (i in 0:n) stock.prices[i+1]<-u^i*(2-u)^(n-i)*S if (type=="euro.call"){ prev<-pmax(stock.prices-K,0) for (k in (n-1):0){ deriv.prices<-matrix(0,k+1,1) for (i in 1:(k+1)) deriv.prices[i]<-(1+r*T/n)^(-1)*((1-q)*prev[i]+q*prev[i+1]) prev<-deriv.prices } } if (type=="euro.put"){ prev<-pmax(K-stock.prices,0) for (k in (n-1):0){ deriv.prices<-matrix(0,k+1,1) for (i in 1:(k+1)) deriv.prices[i]<-(1+r*T/n)^(-1)*((1-q)*prev[i]+q*prev[i+1]) prev<-deriv.prices } } if (type=="amer.call"){ prev<-pmax(stock.prices-K,0) for (k in (n-1):0){ deriv.prices<-matrix(0,k+1,1) stock.prices<-matrix(0,k+1,1) for (i in 0:k) stock.prices[i+1]<-u^i*(2-u)^(k-i)*S for (i in 1:(k+1)) deriv.prices[i]<- max(stock.prices[i]-K,(1+r*T/n)^(-1)*((1-q)*prev[i]+q*prev[i+1])) prev<-deriv.prices } } if (type=="amer.put"){ prev<-pmax(stock.prices-K,0) for (k in (n-1):0){ deriv.prices<-matrix(0,k+1,1) stock.prices<-matrix(0,k+1,1) for (i in 0:k) stock.prices[i+1]<-u^i*(2-u)^(k-i)*S for (i in 1:(k+1)) deriv.prices[i]<- max(K-stock.prices[i],(1+r*T/n)^(-1)*((1-q)*prev[i]+q*prev[i+1])) prev<-deriv.prices } } return(prev[1]) } create.tree<-function(cnum,lnum, p=rep(1/cnum,cnum),rets=c(0.9,1.0,1.1),s0=100) { # cnum = number of childs (2,3,...) # lnum = number of levels (2,3,...) # p is cnum-vector: probability distributio for childs # result: list of # child, sibling, value, begs (lnum-vector of pointers to the begs of levels nodenum<-1 prenum<-1 lnums<-matrix(1,lnum,1) for (i in 1:(lnum-1)){ nodenum<-nodenum+prenum*cnum lnums[i+1]<-prenum*cnum prenum<-prenum*cnum } child<-matrix(0,nodenum,1) sibling<-matrix(0,nodenum,1) value<-matrix(0,nodenum,1) begs<-matrix(0,lnum,) begs[1]<-1 value[1]<-s0 child[1]<-2 sibling[1]<-0 node<-0 for (i in 1:(lnum-1)){ for (j in 1:lnums[i]){ node<-node+1 child[node]<-node+(j-1)*cnum+j if (j/cnum==round(j/cnum)) sibling[node]<-0 else sibling[node]<-node+1 } } tree<-list(child,sibling,value,begs) return(tree) } hedging.data<-function(S,K,t=0.25,r=0, dendat=NULL,periods=1,daysinyear=250,B=NULL,seed=1,method=2,IV=0.2,h=0) { # S is n-vector of stock prices, the last obs is made today # dendat is n times d matrix of explanatory variables # K is the strike price # t is the time to expiration (fractions of year) # r intrest rate (annual) n<-length(S) # dim(dendat)[1] dT<-round(t*daysinyear) # time in days newn<-n-dT S0<-S[n] if (!is.null(dendat)){ d<-dim(dendat)[2] now<-dendat[n,] arg<-matrix(now,d,1) if (d==1) X<-matrix(dendat[1:newn],newn,1) else X<-dendat[1:newn,] } else{ arg<-NULL X<-NULL } if (is.null(B)){ Z<-matrix(0,newn,periods) payout<-matrix(0,newn,1) if (method==1){ for (i in 1:newn){ payout[i,1]<-max(0,S[i+dT]/S0-K/S0) for (jj in 1:periods){ Tbeg<-i Tend<-Tbeg+floor(dT/periods) Z[i,jj]<-S[Tend]/S[Tbeg]-exp(r*t/periods) } } } else{ for (i in 1:newn){ Kadj<-S[i]*K/S0 P<-bs(S[i],Kadj,IV,r,t,type="call") payout[i,1]<-(S0/S[i])*max(0,S[i+dT]-Kadj)#-exp(r*t)*P for (jj in 1:periods){ Tbeg<-i Tend<-Tbeg+floor(dT/periods) Z[i,jj]<-S0*(S[Tend]/S[Tbeg]-exp(r*t/periods)) } } } } else{ R<-returns(S,method="loga") set.seed(seed) booti<-matrix(0,dT,B) for (b in 1:B){ for (i in 1:dT){ res<-ceiling(runif(1)*length(R)) booti[i,b]<-R[res]+h*rnorm(1) } } Z<-matrix(0,B,periods) payout<-matrix(0,B,1) for (i in 1:B){ suma<-sum(booti[,i]) payout[i,1]<-S0*max(0,exp(suma)-K/S0) for (jj in 1:periods){ Tbeg<-1 Tend<-Tbeg+floor(dT/periods)-1 zuma<-sum(booti[Tbeg:Tend,i]) Z[i,jj]<-S0*(exp(zuma)-exp(r*t/periods)) } } } return(list(arg=arg,X=X,Z=Z,payout=payout)) } hedging<-function(payout,Z,parlo=0,parup=1,crit="var",alpha=0.05) { # Z is n times p matrix of reponse variables; linear combination # of the rows of Z will be the final response variable # we find b for the next period n<-dim(Z)[1] p<-dim(Z)[2] if (crit=="var"){ fn<-function(b) { y<-b*Z[,1]-payout return(sum(y^2)/n-(sum(y)/n)^2) } } else if (crit=="quantile"){ fn<-function(b) { y<-b*Z[,1]-payout return(-quantile(y,probs=alpha)) } } else if (crit=="l1"){ fn<-function(b) { y<-b*Z[,1]-payout ka<-sum(y)/n return(sum(abs(y-ka))) } } par<-rep(0.5,p)/p # initial value par.lower<-rep(parlo,p) par.upper<-rep(parup,p) control<-donlp2.control(silent=TRUE) curp<-donlp2(par=par,fn=fn, par.lower=par.lower, par.upper=par.upper, # A=A, lin.upper=lin.upper, lin.lower=lin.lower, control=control) curppi<-curp\$par return(curppi) } markowitz<-function(mu,sigma,gamma,unconstrained=TRUE,lower=0,upper=1,optim=FALSE) { N<-length(mu) I<-diag(1,N) j<-matrix(1,N,1) if (unconstrained){ invsigma<-solve(sigma,I) lambda<-(1-t(j)%*%invsigma%*%mu/gamma)/(t(j)%*%invsigma%*%j/gamma) b<-(1/gamma)*invsigma%*%(mu+as.numeric(lambda)*j) } else{ # constrained fn<-function(b0){ w<-matrix(c(b0,1-sum(b0)),N,1) resu<-t(w)%*%mu-(gamma/2)*t(w)%*%sigma%*%w return(resu) } par<-matrix(rep(1/N,N-1),N-1,1) k<-2*(N-1)+2 # constraint is ui %*% theta - ci >= 0 ui<-matrix(0,k,N-1) ui[1:(N-1),]<-diag(1,N-1) ui[N:(2*(N-1)),]<-diag(-1,N-1) ui[k-1,]<-1 ui[k,]<--1 ci<-matrix(lower,k,1) ci[N:(2*(N-1)),]<--upper ci[k-1,1]<-lower ci[k,1]<--upper if (N>2){ bb<-constrOptim(par,fn,grad=NULL,ui=ui,ci=ci) } else{ # N=2 if (optim){ bb<-optim(par,fn,gr=NULL,method="L-BFGS-B",lower=lower,upper=upper) b<-c(bb\$par,1-sum(bb\$par)) } else{ sigma12<-sigma[1,2] sigma1<-sigma[1,1] sigma2<-sigma[2,2] bb<-gamma^(-1)*(mu[2]-mu[1]-gamma*(sigma12-sigma1))/(sigma1+sigma2-2*sigma12) b<-matrix(0,2,1) b[2]<-min(1,max(0,bb)) b[1]<-1-b[2] } } # end N=2 } # end constrined return(b) } mc.1step<-function(mc,s0=100,r=0,dt=1/250,type="call",moneyness=1) { # mc is T*M matrix of M Monte Carlo samples of length T # moneyness = K/s0 T<-dim(mc)[1] # T=1 M<-dim(mc)[2] dx1<-mc[1,]-s0*exp(r*dt) if (type=="call") H1<-s0*pmax(mc[1,]/s0-moneyness,0) #xi<-cov(c(mc),H1)/var(c(mc)) xi<-cov(dx1,H1)/var(dx1) H0<-exp(-r*dt)*(mean(H1)-xi*mean(dx1)) return(list(H0=H0,xi=xi)) } mc.dyna.2step<-function(mc,s0=100,h=1,r=0,dt=1/250,type="call",moneyness=1) { # mc is T*M matrix of M Monte Carlo samples of length T # moneyness = K/s0 T<-dim(mc)[1] # T<-2 M<-dim(mc)[2] dx2<-mc[2,]-mc[1,]*exp(r*dt) dx1<-(mc[1,]-s0*exp(r*dt))*exp(r*dt) if (type=="call") H2<-s0*pmax(mc[2,]/s0-moneyness,0) edx2<-matrix(0,M,1) edx2sq<-matrix(0,M,1) eH2<-matrix(0,M,1) edx2H2<-matrix(0,M,1) eH2mw1dx2<-matrix(0,M,1) for (i in 1:M){ arg<-mc[1,i] x<-mc[1,] y<-dx2 edx2[i]<-kernesti.regr(arg,x,y,h) edx2sq[i]<-kernesti.regr(arg,x,y^2,h) eH2[i]<-kernesti.regr(arg,x,H2,h) edx2H2[i]<-kernesti.regr(arg,x,dx2*H2,h) } k1<-1-edx2^2/edx2sq H1<-(1/k1)*(eH2-edx2*edx2H2/edx2sq) k0<-mean(k1)-mean(k1*dx1)^2/mean(k1*dx1^2) H0<-(1/k0)*(mean(k1*H1)-mean(k1*dx1)*mean(k1*dx1*H1)/mean(k1*dx1^2)) xi1<-mean((H1-H0)*k1*dx1)/mean(k1*dx1^2) return(list(H0=H0,xi=xi1)) } mc.dyna<-function(mc,s0=100,h=1,r=0,dt=1/250,type="call",moneyness=1) { # mc is T*M matrix of M Monte Carlo samples of length T # moneyness = K/s0 T<-dim(mc)[1] M<-dim(mc)[2] if (type=="call") HT<-s0*pmax(mc[T,]/s0-moneyness,0) ht<-HT kt<-rep(1,M) for (t in (T-1):1){ dx<-mc[t+1,]*exp(r*(t+1)*dt)-mc[t,]*exp(r*t*dt) ek<-matrix(0,M,1) ekdx<-matrix(0,M,1) ekdxsq<-matrix(0,M,1) ekh<-matrix(0,M,1) ekdxh<-matrix(0,M,1) for (i in 1:M){ arg<-mc[t,i] x<-mc[t,] ek[i]<-kernesti.regr(arg,x,kt,h) ekdx[i]<-kernesti.regr(arg,x,kt*dx,h) ekdxsq[i]<-kernesti.regr(arg,x,kt*dx^2,h) ekh[i]<-kernesti.regr(arg,x,kt*ht,h) ekdxh[i]<-kernesti.regr(arg,x,kt*dx*ht,h) } kt<-ek-ekdx^2/ekdxsq ht<-(1/kt)*(ekh-ekdx*ekdxh/ekdxsq) } dx1<-mc[1,]*exp(r*dt)-s0 k0<-mean(kt)-mean(kt*dx1)^2/mean(kt*dx1^2) h0<-(1/k0)*(mean(kt*ht)-mean(kt*dx1)*mean(kt*dx1*ht)/mean(kt*dx1^2)) xi<-mean((ht-h0)*kt*dx1)/mean(kt*dx1^2) return(list(H0=h0,xi=xi)) } mc.hedge<-function(S,K,v,r,t,type="call",lkm=1000, mu=0,sig=1,annu.mu=NULL, dist="norma",df=4,model.type="incre", seed=0,daysinyear=250,n=round(t*daysinyear),freq=round(daysinyear*t/n),steps=1) { payout<-matrix(0,lkm,1) if (steps==1){ X1<-matrix(0,lkm,1) for (ii in 1:lkm){ seedii<-seed+ii ts<-sim.ts(n,S0=S, mu=mu,sig=sig,annu.mu=annu.mu,annu.vola=v, dist=dist,df=df,model.type=model.type,seed=seedii, freq=freq,daysinyear=daysinyear) S1<-ts[length(ts)] X1[ii]<-S0*(S1/S0-exp(r*t)) payout[ii]<-max(0,S1-K) } delta<-cov(payout,X1)/var(X1) price<-exp(-r*t)*mean(payout-delta[1,1]*X1) delta2<-0 } if (steps==2){ X1<-matrix(0,lkm,1) X2<-matrix(0,lkm,1) for (ii in 1:lkm){ seedii<-seed+ii ts<-sim.ts(n,S0=S, mu=mu,sig=sig,annu.mu=annu.mu,annu.vola=v, dist=dist,df=df,model.type=model.type,seed=seedii, freq=freq,daysinyear=daysinyear) S2<-ts[length(ts)] S1<-ts[round(length(ts)/2)] X1[ii]<-S0*(S1/S0-exp(r*t)) X2[ii]<-S1*(S2/S1-exp(r*t)) payout[ii]<-max(0,S2-K) } delta<-cov(payout,X1)/var(X1) delta2<-cov(payout,X2)/var(X2) price<-exp(-r*t)*mean(payout-delta[1,1]*X1-delta2[1,1]*X2) } return(list(price=price,hedge=delta,hedge2=delta2)) } mc.local.2step<-function(mc,s0=100,h=1,r=0,dt=1/250,type="call",moneyness=1) { # mc is T*M matrix of M Monte Carlo samples of length T # moneyness = K/s0 T<-dim(mc)[1] # T<-2 M<-dim(mc)[2] S1<-mc[1,] S2<-mc[2,] if (type=="call") H2<-s0*pmax(mc[2,]/s0-moneyness,0) eh2<-matrix(0,M,1) es2<-matrix(0,M,1) es2sq<-matrix(0,M,1) es2h2<-matrix(0,M,1) for (i in 1:M){ arg<-mc[1,i] x<-mc[1,] eh2[i]<-kernesti.regr(arg,x,H2,h) es2[i]<-kernesti.regr(arg,x,S2,h) es2sq[i]<-kernesti.regr(arg,x,S2^2,h) es2h2[i]<-kernesti.regr(arg,x,S2*H2,h) } covs2h2<-es2h2-es2*eh2 vars2<-es2sq-es2^2 b2<-covs2h2/vars2 a2<-eh2-b2*es2 #H1<-exp(-r*dt)*eH2+b2*(S1-exp(-r*dt)*e1s2 H1<-exp(-r*dt)*a2+b2*S1 b1<-cov(S1,H1)/var(S1) a1<-mean(H1)-b1*mean(S1) #H0<-exp(-r*dt)+mean(H1)+b1*(s0-exp(-r*dt)*mean(S1)) H0<-exp(-r*dt)*a1+b1*s0 xi1<-b1 return(list(H0=H0,xi=xi1)) } sharpe.confidence<-function(x,alpha=0.05,ker=sqrt(250)) { mu<-mean(x) sigma<-sd(x) sharpe<-mu/sigma T<-length(x) ST<-c(mu,mean(x^2)) x1<-ST[1] x2<-ST[2] grad<-c(x2/(x2-x1^2)^(3/2),-0.5*x1/(x2-x1^2)^(3/2)) Z<-matrix(0,T,2) Z[,1]<-x Z[,2]<-x^2 Psi0<-t(Z-ST)%*%(Z-ST)/T hats<-sqrt(t(grad)%*%Psi0%*%grad) U<-sharpe+T^(-1/2)*qnorm(1-alpha/2)*hats L<-sharpe-T^(-1/2)*qnorm(1-alpha/2)*hats U<-ker*U L<-ker*L sharpe<-ker*sharpe return(list(sharpe=sharpe,L=L,U=U)) } sim.price<-function(retu,s0=100,type="log") { n<-dim(retu)[1] M<-dim(retu)[2] S<-matrix(0,n,M) if (type=="log"){ S[1,]<-s0*exp(retu[1,]) if (n>=2) for (i in 2:n) S[i,]<-S[i-1,]*exp(retu[i,]) } if (type=="gross"){ S[1,]<-s0*retu[1,] if (n>=2) for (i in 2:n) S[i,]<-S[i-1,]*retu[i,] } if (type=="net"){ S[1,]<-s0*(retu[1,]+1) if (n>=2) for (i in 2:n) S[i,]<-S[i-1,]*(retu[i,]+1) } return(S) } sim.return<-function(n,M,type="gauss",mu=0,sig=1,seed=1) { set.seed(seed) if (type=="gauss") retu<-matrix(mu+sig*rnorm(n*M),n,M) return(retu) } tail.hill<-function(data,p=0.1,m=round(p*length(data)),type="left",h=NULL) { if (type=="left") or<-order(data) else or<-order(data,decreasing=TRUE) if (is.null(h)){ dat<-data[or] dit<-dat[1:m] c<-dit[m] alpha<-m/sum(log(dit/c)) } else{ n<-length(data) x<-seq(1,n) arg<-n w<-kernesti.weights(arg,x,h=h) we<-w[or] cs<-cumsum(we) ind<-(cs<=p) if (sum(ind)==0) m<-1 else m<-max(x[ind]) wet<-we[1:m] dat<-data[or] dit<-dat[1:m] c<-dit[m] alpha<-p/sum(wet*log(dit/c)) } return(alpha) } tail.index<-function(data,p=0.1,m=round(p*length(data)),type="left",h=NULL) { n<-length(data) if (is.null(h)){ if (type=="left"){ or<-order(data) dat<-data[or] dit<-dat[1:m] x<--log(-dit) y<-log(seq(1,m)/n) } else{ or<-order(data,decreasing=TRUE) dat<-data[or] dit<-dat[1:m] x<--log(dit) y<-log(seq(1,m)/n) } alpha<-cov(x,y,use="complete.obs")/var(x,na.rm=TRUE) } else{ if (type=="left") or<-order(data) else or<-order(data,decreasing=TRUE) xs<-seq(1,n) arg<-n w<-kernesti.weights(arg,xs,h=h) we<-w[or] cs<-cumsum(we) ind<-(cs<=p) if (sum(ind)==0) m<-1 else m<-max(xs[ind]) wet<-we[1:m] dat<-data[or] dit<-dat[1:m] if (type=="left") x<--log(-dit) else x<--log(dit) y<-log(seq(1,m)/n) mx<-sum(wet*x,na.rm=TRUE)/sum(wet,na.rm=TRUE) my<-sum(wet*y,na.rm=TRUE)/sum(wet,na.rm=TRUE) mxy<-sum(wet*x*y,na.rm=TRUE)/sum(wet,na.rm=TRUE) mx2<-sum(wet*x^2,na.rm=TRUE)/sum(wet,na.rm=TRUE) alpha<-(mxy-my*mx)/(mx2-mx^2) } return(alpha) }