@*****************************************************@ @ FFT.gss - supplementary program to @ @ @ @ "Introduction to FFT in Finance" @ @ @ @ (c) 2004 Ales Cerny, all rights reserved @ @*****************************************************@ /* This program implements binomial option pricing model using fast Fourier transform */ /* Comments on commands in this program hsec time from midnight in hundreds of a second seqa(start,incr,no) arithmetic progression with starting point start, increment incr and number of elements no seqm(start,quot,no) multiplicative (geometric) progression with quotient quot ceil(X) the nearest bigger integer (ceiling) to X X' matrix X transposed X.*Y matrix X (m x n) multiplied 'element by element' by matrix Y Y must be either (m x n) or (m x 1) or (1 x n) the result is always (m x n) 0|1|2 a column vector with elements 0,1,2 0~1~2 a row vector with elements 0,1,2 zeros(m,n) (m x n) matrix of zeros ones(m,n) (m x n) matrix of ones rows(X) row dimension of matrix X cols(X) column dimension of matrix X minc(X) returns the smallest value in each column of X as a column vector, i.e. maxc((1~2~3)|(4~0~0)) = 1|0|0 maxc(X) returns the largest value in each column of X as a column vector, i.e. maxc((1~2~3)|(4~0~0)) = 4|2|3 for tt (st,end,incr) loop over tt starting at st, increasing tt each time by incr, stop endfor when tt is greater than end */ new; @ clear the workspace @ cls; @ clear screen @ @***************************@ @ trading time @ @ parameters @ @***************************@ Minute = 1; Hour = 60; HoursInDay = 8; DaysInWeek = 5; DaysInMonth = 21; Day = HoursInDay*Hour; Week = DaysInWeek*Day; Month = DaysInMonth*Day; Year = 12*Month; @***************************@ @ Hedging Parameters @ @***************************@ T=3*month; @ Time to maturity @ RehedgeInterval=minute/6; @ Trading period @ S0=5100; @ Initial stock price @ strike=5355; @***************************@ @ Transformation of @ @ log returns @ @***************************@ UnitTime=month; R1safe=1.0033; @ monthly safe return @ R1=1.053~0.965; @ monthly return @ PDistr=0.5~0.5; @ prob. density of monthly returns @ lnR1=ln(R1); @ monthly log return @ mu1=lnR1*PDistr'; @ expected monthly log return @ sig1=sqrt(((lnR1-mu1)^2)*PDistr'); @ volatility of monthly log return @ dt = RehedgeInterval/UnitTime; lnRdt=mu1*dt+(lnR1-mu1)*sqrt(dt); @ log return over rehedging interval @ Rdt=exp(lnRdt); Rdtsafe=R1safe^dt; @************************@ @ Risk-neutral @ @ probabilities @ @************************@ QDistr=((Rdtsafe-Rdt[2])~(Rdt[1]-Rdtsafe))/(Rdt[1]-Rdt[2]); @************************@ @ grid indexation @ @************************@ Tidx=ceil(T/RehedgeInterval)+1; @ Number of trading dates @ dlnS=lnRdt[1]-lnRdt[2]; @ increment on log price grid @ highlnRdt=lnRdt[1]; @ the highest return over one period @ @ there are tt live cells at time tt, highest stock price at the top @ @ log price at cell 1 at time tt is ln(S0)+(tt-1)*highlnRdt @ @ log price at cell ii at time tt is ln(S0)+(tt-1)*highlnRdt-(ii-1)*dlnS @ @************************@ @ option payoff @ @************************@ S_T=seqa(ln(S0)+(Tidx-1)*highlnRdt,-dlnS,Tidx); @ log price at maturity @ S_T=minc(exp(S_T)'|(ones(1,Tidx)*100*S0)); @ stock price at maturity @ C=maxc((S_T-strike)'|zeros(1,rows(S_T))); @ option payoff at maturity @ @************************@ @ pricing @ @************************@ starttime = hsec; @ start of computation @ pkernel=QDistr'/Rdtsafe|zeros(Tidx-2,1); @ pricing kernel @ CharFction=(nextn(Tidx)*fftn(pkernel)); Cimage=ffti(C); if RehedgeInterval>1*minute; for ii (1,100,1); C = fftn( Cimage.*(CharFction^(Tidx-1)) ); Cnext = fftn( Cimage.*(CharFction^(Tidx-2)) ); endfor; runtime=(hsec - starttime)/10000; else; C = fftn( Cimage.*(CharFction^(Tidx-1)) ); Cnext = fftn( Cimage.*(CharFction^(Tidx-2)) ); runtime=(hsec - starttime)/100; endif; print "______________________________________________"; print " "; format /rd 12,2; print "trading frequency in minutes " RehedgeInterval; print " "; format /rd 12,8; print "binomial FFT price " real(C[1]); print "Black-Scholes price " EuropeanBSCall(S0,strike,ln(R1safe),0,T/UnitTime,sig1); print " "; format /rd 12,8; print "binomial delta " real((Cnext[1]-Cnext[2])/S0/(Rdt[1]-Rdt[2])); print "Black-Scholes delta " cdfn((ln(S0/strike)+(sig1^2/2+ln(R1safe))*T/UnitTime)/(sig1*sqrt(T/UnitTime))); print " "; format /rd 12,3; print "required time in seconds " runtime; print "______________________________________________"; format /rd 16,8;