Chapter 1 基本…下行继续区分大小写变量合法命名：只含字母、数字、下划线，并以字母开头who显示定义变量whos显示定义变量及其信息clear all 清除所有存在变量help “?”查询函数和运用方法lookfor “?”查询与关键词有关的所有函数Integer整数real实数complex复数inf无限大NaN非数字format long显示多位小数short 4位bank 2位round（）四舍五入fix()去小数部分floor()下取ceil()上取基本运算+ - * / ^ pi sqrt（x）开方exp（x）e的x次方prod() 连乘abs（）求绝对值max（）最大值min（）最小值log（）自然对数log2（）log10()sin() cos() tan() cot() asin() acos() atan()矩阵、向量矩阵中空格或“，”分隔列“；”或回车分隔行transpose()或矩阵’矩阵转置inv（）矩阵的倒置diag（）将某矩阵除对角线元素外，其他元素为0，形成一个新的矩阵ones（2,3）两行三列的全是1的矩阵zeros（3,4）三行四列的全是0的矩阵eye（4）4行4列单位矩阵q13(2,4)矩阵q13第二行第四列数据q13(:,3) 矩阵q13第3列数据q13(2,:) 矩阵q13第2行数据q13(2:3,2:3) 矩阵q13 2到3行，2到3列的数据+ - 同行列数的矩阵才能加减*一般矩阵乘法.*对应元素相乘同理./ .^ a：b：c 以a开始，以b为间隔，最大数小于c 的数列字符‘’括起来A为矩阵sum（A）对列求和（得行）sum（A，2）对行求和（得列）sum(sum（A）)对列求和再对行求和（得一个数）prod（A）对列连乘min（A）输出每一列的最小值B为向量[m,j]=max（B）m=最大值取值，j=最大值位置[m,j]=min（B）m=最小值取值，j=最小值位置逻辑运算判断：a==b a等于ba~=b a不等于b < > <= >=正确输出1，错误输出0画图plot（x,y）plot(x,y, '（颜色）（形状）') linspace（a,b,n）把起点为a终点为b的直线等分为n份title('Graph of y = sin(3\pix)') title图表名称xlabel('x axis') x轴名称ylabel('sin(3\pix)') y轴名称plot(x,y,'.',x,cos(3*pi*x),'g*')两条线画于一图legend('Sin curve', 'Cos curve')添加图标grid 加格子hold on 保留原曲线，可用于画多线在一图subplot(121), plot(x,y) 将作图区域分为1行2列，作第1个区域的图相关程序Plot a graph of y = sin(3px) for 0 ≤x ≤1 >> N = 100; h = 1/N;>> x = 0:h:1;>> y = cos(3*pi*x);>> plot(x,y)The command linspace (a,b,n) generates n + 1 equispaced points between a and b>> x = linspace (0,1,101);>> y = sin(3*pi*x);>> plot(x,y)随机数rand 0到1的随机数rand（2,3）2行3列随机矩阵rand（3）3行正矩阵随机矩阵randn 标准正态分布随机数y = random('name',A1,A2,A3,m,n)r5 = random('Normal',0,1,2,4)r6 =random('Poisson',6,1,6)文档读写IBM = xlsread('IBM.xls');IBM = textread('IBM.txt');IBM = dlmread('IBM.txt','',1,1);IBM = importdata('IBM.xls');m文件which mfile（文件名）查看路径edit mfile 编辑脚本mfile与函数mfile的区别，前者是全局变量，后者是局部变量function [var1 var2 …] =functionname(arg1,arg2, …)循环for…end while…end if…end if…elseif…elseif …else…end for…for…end…endswitch case… case… otherwise… end相关程序For Loops画多个图形>> x = -1:.05:1;>> for n = 1:8ubplot(4,2,n);plot(x,sin(n*pi*x));end求和>> sum = 0>> for i = 1:100sum = sum + i;end>> sumWhile LoopsWhat is the greatest value of n that can be used in the sum1^2 + 2^2 + … + n^2and get a value of less than 100?>> S = 1;>> n = 1;>> while S+ (n+1)^2 < 100n = n+1;S = S + n^2;end>> [n, S]If...Then Loop>> S = 1; n = 1;>>for i=1:100if S+(n+1)^2 < 100S = S+ (n+1)^2; n = n+1;endend>> [n, S]Chapter 2 现金流分析pvvar求变动现金流量的现值pvfix求固定现金流量的现值fvvar求变动现金流量的终值fvfix求固定现金流量的终值PV = pvvar(CashFlow, Rate, CFDates)程序>> CashFlow = [100 300 450]>> CFDates = ['01/12/1987';'02/14/1988';'03/03/1988']>> PV = pvvar(CashFlow, 0.09, CFDates)注意日期是列数据>> pvvar([-15000 3000 4500 5000 6800], 0.08) ans = 603.1667没有现金流日期时，默认为没两个现金流相隔一年Irr固定周期的内部收益率Xirr变动周期的内部收益率Return = irr(CashFlow)Return = xirr(CashFlow, Dates)程序>> CF = [-10000 3000 4500 5000];>> Return = irr(CF)Return =0.1105>> Dates=['01/12/00'; '02/14/01'; '09/03/01';'12/31/02'];>> Return = xirr(CF, Dates)Return =0.1177cfdur and cfconv 久期和凸性[Dur, ModDur] = cfdur(CashFlow, Yield) Conv = cfconv(CashFlow, Yield)程序>> CashFlow=[5 5 5 5 5 105];>> [Dur, ModDur] = cfdur(CashFlow, 0.05)Dur =5.3295ModDur =5.0757>> Conv = cfconv(CF, 0.05)Conv =95.5410债券价格和收益率[P, I应付利息] = bndprice(Yield, CouponRate票面利率, Settle交收日, Maturity到期) 日期一定要加单引号p=price+accruedintYield = bndyield(Price, CouponRate, Settle, Maturity)Duration = bnddurp(Price, CouponRate, Settle, Maturity)Duration = bnddury(Yield, CouponRate, Settle, Maturity)Convexity = bndconvp(Price, CouponRate, Settle, Maturity)Convxity = bndconvy(Yield, CouponRate, Settle, Maturity)程序Price-yield curve>> yields = 0.01: 0.01: 0.20;>> [P I] = bndprice(yields, 0.1, '08/10/07','12/31/20');>> plot (yields, P+I);>> grid on;>> xlabel('Yield');>> ylabel('Price') ;>> Title('Price-Yield Curve'); 利率免疫duration=holding period% SET BOND FEATURES (bondimmun.m) >> settle = '28-Aug-2007';>> maturity = ['15-Jun-2012' ; '31-Oct-2017' ;'01-Mar-2027'];>> couponRate = [0.07 ; 0.06 ; 0.08];>> yield = [0.06 ; 0.07 ; 0.075];% COMPUTE DURATIONS AND CONVEXITIES>> duration = bnddury(yield, couponRate, settle, maturity);>> convexity = bndconvy(yield, couponRate, settle, maturity);% COMPUTE PORTFOLIO WEIGHTS>> A = [duration'; convexity'; 1 1 1];>> b = [ 10; 160; 1];>> weights = A\b有效前沿[PortRisk, PortReturn, PortWts] =frontcon(ExpReturn历史收益率, ExpCovariance历史协方差, NumPorts返回结果的个数, PortReturn目标收益, AssetBounds资产界限, Groups, GroupBounds组合界限)程序>> r = [0.2 0.1];两种资产的历史收益率>> s = [0.2 -0.1; -0.1, 0.4];两种资产的历史协方差>> [Risk, Return, Wts] = frontcon(r, s, 5);>> [Risk, Return, Wts]ans =0.2958 0.1625 0.6250 0.37500.3075 0.1719 0.7187 0.28130.3400 0.1812 0.8125 0.18750.3883 0.1906 0.9062 0.09380.4472 0.2000 1.0000 0最后两列为权重含约束条件的有效前沿[PortRisk, PortReturn, PortWts] =portopt(Return, Cov, [], PortReturn, ConSet)约束条件ConSet = portcons('Default', Num, 'AssetLims', Min, Max,'GroupLims', Group, GroupMin, GroupMax);程序>> r = [0.1 0.2 0.15];>> s = [0.0100 -0.0061 0.0042; -0.0061 0.0400 -0.0252; 0.0042 -0.0252 0.0225];>> ObjReturn = 0.17;>> [Risk, Return, Wts] = portopt(r, s, [], ObjReturn)Return =0.1700Wts =0.0278 0.4278 0.5445>> Return = [0.1 0.2 0.15];>> Cov = 0.01*[0.5 -1 0.4; -1 4 -0.2; 0.4 -0.2 2.3];>> PortReturn = [0.15 0.16];目标收益率在0.15~0.16之间>> Num = 3;>> Min = [0.20 NaN NaN]; 投资组合中每种资产的下限>> Max = [0.5 0.5 0.5]; ……上限>> Group = [1 1 0; 0 1 1];分为两组，第一组由第1、2种资产构成；第二组由……>> GroupMin = [0.2, 0.2]; 第一组的投资比例下限为0.2；第二组……>> GroupMax = [0.8, 0.8]; 第一组的投资比例上限为0.8；第二组……>> ConSet = portcons('Default', Num,'AssetLims', Min, Max,'GroupLims', Group,GroupMin, GroupMax);>> [PortRisk, PortReturn, PortWts] =portopt(Return, Cov, [], PortReturn, ConSet) PortRisk = 0.0724 0.0910 PortReturn = 0.1500 0.1600 PortWts =0.4000 0.4000 0.2000 0.2606 0.4606 0.2789可借贷无风险资产时的有效前沿 [RiskyRisk, RiskyReturn, RiskyWts,RiskyFraction 所有风险资产占总资产的比例, OverallRisk, OverallReturn] =portalloc (PortRisk, PortReturn, PortWts,RisklessRate, BorrowRate, RiskAversion 风险厌恶程度) 程序>> Ret = [0.1 0.2 0.15];>> Cov = [ 0.005 -0.010 0.004; -0.010 0.040 -0.002; 0.004 -0.002 0.023];>> [Risk, Return, Wts] = frontcon(Ret, Cov, 20); >> RiskfreeR = 0.08; BorrowR = 0.12; >> RiskAversion = 3;>> portalloc (Risk, Return, Wts, RiskfreeR, BorrowR, RiskAversion); (return the graphic)>>[RiskyRisk, RiskyReturn, RiskyWts,RiskyFraction, OverallRisk, OverallReturn] = portalloc (Risk, Return, Wts, RiskfreeR, BorrowR, RiskAversion)最优化问题[x, fval] = fmincon (fun 目标函数, x0从x0开始试值, A, b, Aeq, beq, lb, ub, nonlcon 非线性约束,options)options = optimset ('LargeScale','off') 注意：x1要加括号x （1）；目标函数要建立m 文件；非线性约束要建立m 文件程序Step 1: Write an M-file Myobj.m.function function [f] = Myobj(x) f =exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+x(2)+1);Step 2: Write an M-file Mycon.m for the constraintsfunction [c, ceq] = Mycon(x)% Nonlinear inequality constraintsc = [1.5 + x(1)*x(2) - x(1) - x(2); -x(1)*x(2) - 10];% Nonlinear equality constraintsceq = [];Step 3: Invoke constrained optimization routineA = [1 1; 2 -1]; b = [20; 10]; Aeq = []; beq = [];lb = [-10 -10]; ub = [10, 10];x0 = [-1,1]; % Make a starting guess at the solutionoptions = optimset('LargeScale','off');[x, fval] = fmincon(@Myobj, x0, A, b, Aeq, beq, lb, ub, @Mycon, options)B-S 模型[Call, Put] = blsprice (Price, Strike, Rate, Time, V olatility, Yield)[CallDelta, PutDelta] = blsdelta (Price, Strike, Rate, Time, V olatility, Yield)Gamma = blsgamma (Price, Strike, Rate, Time, V olatility, Yield)[CallEl, PutEl] = blslambda (Price, Strike, Rate, Time, V olatility, Yield)[CallRho, PutRho]= blsrho (Price, Strike, Rate, Time, V olatility, Yield)[CallTheta, PutTheta] = blstheta (Price, Strike, Rate, Time, V olatility, Yield)Vega = blsvega (Price, Strike, Rate, Time, V olatility, Yield)Volatility = blsimpv (Price, Strike, Rate, Time, Value, Limit, Yield, Tolerance, Class)Chapter 3 二叉树模型European and American CallEuropean PutAmerican Put 程序欧式二叉树function [] = BinomialEuro()T=4/12 % experiation daten=40; %total number of periods K=300; %exercise pricer=0.08;%risk free interest rate sigma=0.3;%vloatility of S0=300; deltat=T/n;u=exp(sigma*sqrt(deltat)); d=exp(-sigma*sqrt(deltat));p=(exp(r*deltat)-d)/(u-d); %risk adjusted probability q=1-p;for i=1:n+1s(i)=u^(n+1-i)*d^(i-1)*S0; % stock price w(i)=nchoosek(n,i-1)*p^(n+1-i)*q^(i-1); % probabilityx(i)=max(s(i)-K,0); % call: min ()..:()0 ()0xobj f x s t c x ceq x A x bAeq x beq lb x ub ≤=⋅≤⋅=≤≤非线性不等式约束 非线性等式约束线性不等式约束 线性等式约束 上下界 1221122121212121212: min (4241)..: 1.50 100 20 210-10,10x x obj ex x x x xs t x x x x x x x x x x x x ++++--+≤--≤+≤-≤≤≤(),r q t e dp u e d e u d -∆--⇒===-0,0,,0max{,0)TrT T i T i i f e p S K -==-∑0,0,,0max{,0)TrT T i T i i f e p K S -==-∑,01,1,1max{max{,0}, ()}t i i r tt i t i t i f K S u d e pf qf --∆+++=-+y(i)=max(-s(i)+K,0); % European putcv(i)=w(i)*x(i);pv(i)=w(i)*y(i);endC=sum(cv(:))*exp(-r*T)P=sum(pv(:))*exp(-r*T)美式二叉树function [c] = BinomialAmer()T=4/12; % experiation daten=4; %total number of periodsK=300; %exercise pricer=0.08;%risk free interest ratesigma=0.3;%vloatility ofS0=300;deltat=T/n;u=exp(sigma*sqrt(deltat));d=exp(-sigma*sqrt(deltat));p=(exp(r*deltat)-d)/(u-d); %risk adjusted probabilityq=1-p;s(1,1)=S0; %first stage price of stock;for t=2:n+1 %stock price from stage 2 to n+1 for i=1:ts(t,i)=u^(t-i)*d^(i-1)*S0;endendfor i=1:n+1 %option price at last stagex(n+1,i)=max(K-s(n+1,i),0);endfor t=n:(-1):1for i=1:tx(t,i)=(q*x(t+1,i)+p*x(t+1,i+1))/exp(r*deltat); % option price at node (t,i) without exercisex(t,i)=max(x(t,i),K-s(t,i)); %check for exerciseendendc=x(1,1);Chapter 4 蒙特卡罗模拟计算面积Input: Total number n of random points Output: AREA = approximate area under the curveStep 1: Initialize初始化: COUNTER计数器= 0 Step 2: For i= 1, 2, …, n, do Steps 3-5Step 3: Generate random coordinates坐标xi and yi, satisfying a≤xi≤b and 0 ≤yi≤M.Step 4: Calculate f (xi) for the random xiStep 5: If yi≤f (xi), increase the COUNTER by 1. Otherwise, leave COUNTER as is.Step 6 CalculateAREA = M(b –a) COUNTER / nStep 7 OUTPUT(AREA), STOP程序%The area under the curve y = cos x over the interval [-pi/2, pi/2]function [area] = MC1(n)X=rand(n,1)*pi-pi/2 %generating n uniformly distributed RVs over the interval [-pi/2, pi/2]Y=rand(n,1) %generating n uniformly distributed RVs over the interval [0,1] Z=cos(X); %calculate the value of f(x)counter=0;for i = 1 : nif Y(i) <= Z(i)counter=counter+1;endendarea=pi*counter/n;估计pi程序function [Mypi] = Mypi()n = 10000; %number of replicatescount = 0;for i = 1:nx = rand;y = rand;if x^2+y^2 <= 1count = count + 1;endendMypi = 4*count/n;计算积分程序%estimate the integral积分of sin(x)cos(x) function [Myinte] = Myinte()n = 10000; %number of replicatesa = 1; % lower boundb = 3; % upper boundcount = 0;for i = 1:nx = 2*rand+1;count = count + cos(x)*sin(x);endMyinte = (b-a)*count/n;欧式看涨期权定价亚式看涨期权定价1()[()]()[()]()b b Niia ab af x dx E f x dx b a E f x f xN=-≈=-≈1()[()]()[()]()b b Niia ab af xdx E f x dx b a E f x f xN=-≈=-≈∑⎰⎰()()()2211,,()(0)()ˆiir TirTi in nfor i ngenerate Zset S T S eset c e S TKset c c c nσ-++-===-=+()()()()2121111,,1,,()()()()()ˆj j ijijj jmrTin nr t tfor i nfor j mgenerate Zset S t S t eset S S t S tset c e S Kset c c c nσ---+--+====++=-=+Chapter 5 随机数生成器线性同余产生器程序function[]=line(a,m,n) x(1)=2 for i=1:nx(i+1)=mod(a*x(i),m); u(i+1)=x(i+1)/m; end x u%余数mod （被除数，除数）混合线性同余产生器程序function [u] = lcg(seed, n)% the equation to generate random number is: x(n+1)=(a*x(n)+b) mod m% seed is used to generate the first random number% n is the number of random number wanteda = 1229;b = 1;m = 32768;x = zeros(1,n);x(1)=mod(a*seed+b, m);for i = 2:nx(i)= mod(a*x(i-1)+b, m); endu = x./m; % make the random unmber between 0 and 1plot(u(1:n-1), u(2:n), '.'); % plot the pair of point (u(n), u(n+1))title('Autocorrelaton of LCG'); 计算循环的周期function [p]=lcg(a,b,m,x0,n) x(1)=mod((a*x0+b),m) u(1)=x(1)/m for i=2:nx(i)=mod((a*x(i-1)+b),m) u(i)=x(i)/m if u(i)==u(1) p=i end end洗牌程序function [] = lcgsh(seed, n)% This function implements a shuffling scheme on the numbers generated by% seed is used to generate the first random number% n is the number of random number wantedu = lcg(seed, n+20);Seq = ceil(20*lcg(10, n))+1; w = zeros(1,n);for i = 1:nh = Seq(i); w(i) = u(h); u(h) = u(i+20); endplot(w(1:n-1), w(2:n), '.'); % plot the pair of point (u(n), u(n+1))title('Autocorrelaton of LCG');%% the simplest way % u = lcg(seed, n);% Seq = ceil((n-1)*lcg(10, n))+1; % w = zeros(1,n); %% for i = 1:n% h = Seq(i); % w(i) = u(h); % end %% plot(w(1:n-1), w(2:n), '.'); % plot the pair of point (u(n), u(n+1))% title('Autocorrelaton of LCG');掷骰子程序function [outcome] = dice(n)% this function is a Random Dice Generator % n is the number of trialsU = rand(1, n);outcome = zeros(1,n);for i = 1:nif (0 <= U(i) & U(i) <= 1/6) outcome(i) = 1;elseif (1/6 < U(i) & U(i) <= 2/6) outcome(i) = 2;elseif (2/6 < U(i) & U(i) <= 3/6) outcome(i) = 3;elseif (3/6 < U(i) & U(i) <= 4/6) outcome(i) = 4;elseif (4/6 < U(i) & U(i) <= 5/6) outcome(i) = 5; elseoutcome(i) = 6; end end生成符合possion 过程的随机数Step 1: Set s 1=0, s 2=p 1, j =0. And generate a uniform random number U ;Step 2: Check If s 1≤U ≤s 2, then x =j . And exit. Step 3: Else update s 1=s 2, s 2=s 2+p j+1, j =j +1; Go to Step 2. 程序function [] = Poisson(seed, lambda, n)% This function implements a Poisson random number generator% seed is used to generate the first random number% lambda is the parameter of Poisson distribution % n is the number of random number wanted%calculate the prob. of p(x) for x = 1:m, m is the assumed maximum of x m = 1000; P = ones(1,m); for i = 1:m P(i) =exp(-lambda)*lambda^(i-1)/factorial(i-1); endpoisson = ones(1,n);111mod /i i i i x ax m u x m +++==111()mod /i i i i x ax b m u x m +++=+=rand('state',seed) for i = 1:n s1 = 0; s2 = P(1); U = rand; for j = 1:mif (s1 <= U & U <= s2) poisson(i) = j-1; break elses1 = s2;s2 = s2 + P(j+1); end end endmean(poisson) var(poisson)生成符合连续概率分布的随机数 Draw a random number U ~U (0,1) Return X = F−1 (U ) Example Exponential CDF:Inverse of the CDF: 程序>> rand(…state',1); >> -2*log(rand) ans =1.34Box-Mueller 方法生成符合正态分布的随机数 程序function [w] = Boxmuller(seed, n)% Box-Muller method to generate normally distributed random numbersu = lcg(seed,n); % u = rand(1,n);G1 = zeros(1, n); G2 = zeros(1, n); w = zeros(1, 2*n);for i = 1: n,G1(i) = sqrt(-2 * log(u(i))) * cos(2 * pi * u(n-i+1));G2(i) = sqrt(-2 * log(u(i))) * sin(2 * pi * u(n-i+1)); endw = [G1 G2];plot(G1, G2, '.');title('Plot of G1, G2'); xlabel('G1'); ylabel('G2');[H,P,KSSTA T,CV] = kstest(w,[],0.05,0)figurenormplot(w);Chapter6 减小方差的方法控制变量法程序function []= ContrVariate()clear allS0=90; % inital price of the underlying K=90; % strike price r=0.05; % risk-free ratesigma=0.3; % volatility of the underlying T=0.25; % experiation daten=10000; % number of replications% compute call option price: uncontrolled estimator for i=1:n% simulate stock price using Continous risk neutral dynamics when deltat=0.1deltaw = random('Normal',0,sqrt(T),1,1); % independent increament ST(i) =S0*exp((r-1/2*sigma^2)*T+sigma*deltaw); % stock price at expiration dateY1(i)=max(ST(i)-K,0)*exp(-r*T); % payoff of European Call end% compute call option price: controlled estimator covM=cov(ST,Y1); % covarance coefficient between ST and YcorrM=corrcoef(ST, Y1); % correlation coefficient between ST and Yb=-covM(1,2)/covM(1,1); % the optimal coefficient for controlled variate estimatorfor i=1:nY2(i)=Y1(i)+b*(ST(i)-S0*exp(r*T)); % payoff of European Call using control variate endEurC1=mean(Y1(:)) % uncontrolled estimator EurC2=mean(Y2(:)) % controlled estimatorVarEurC1=var(Y1); % variance of uncontrolled estimatorVarEurC2=var(Y2); % variance of controlled estimatorratio=VarEurC2/VarEurC1 % the ration of the variance of the optimally controlled estimator to the uncontrolled estimator is对偶变量法Step 1: generate random number e 1~U(0,1), andcompute f 1 = g (e 1)Step 2: let e 2 =1-e 1, and compute f 2 = g (e 2)Step 3: compute Example>> rand('seed',0)>> U=rand(1,100);>> X=exp(U); >> I=mean(X) I =1.7461>> S2=var(X) S2 = 0.2499Chapter7 布朗运动/()1, 0x P F x e x θ-==->log(), (0,1)X P P U θ=- *(,)()Cov X Y b Var X =-()([])i i i Y b Y b X E X =--11var()var()2f f <122f f f +=11 1.7183x I e dx e ==-≈⎰>>rand('seed',0) >> U1=rand(1,50); >> U2=1-U1; >> X1=exp(U1); >> X2=exp(U2);>> X=0.5*(X1+X2); >> I=mean(X)I >> =1.7175>> S2=var(X)S2 = 0.0044程序%Generating Brownian Motion Sample Paths function BrownianMotion()n=100000; deltat=1/n;mu=0.1; %mean per yearsigma=0.2; %standard deviation per year X=zeros(4,n);randnumber=normrnd(0,1,4,n);%正态分布（均值，标准差，4行，n 列）for i=1:nX(1,i+1)=X(1,i)+sqrt(deltat)*randnumber(1,i); % standard BWX(2,i+1)=X(2,i)+mu*deltat+sigma*sqrt(deltat)*r andnumber(2,i); % general BWX(3,i+1)=X(3,i)+mu*deltat+sigma*X(3,i)*sqrt(d eltat)*randnumber(3,i); % Ito process with drift=mu, diffussion=sigma*XX(4,i+1)=X(4,i)+mu*X(3,i)*deltat+sigma*X(4,i)*sqrt(deltat)*randnumber(4,i); % Ito process with drift=mu*X, diffussion=sigma*X end跳跃过程Step 1: generate Z ~N (0,1);Step 2: generate N ~Poisson(lD t ); if N =0, set M =1 and go to step 4Step 3: generate Y 1,…Y N from their common distribution and set M =Y 1×…×Y Step4: set程序function [c] = jumpdiffusion()n=10; %total number of trading days r=0.1; %rate of return per yearsigma=0.4612; %standard deviation of the underlyingS0=90.44; %price of underlying;lambda=1*4; %Possion process with lambda per yeardeltat = 1/252; %delta t per step T=10/252; %maturity timew = 1000; %# of replications for s = 1:wX = random('Poisson',T*lambda);%compute the number of jumps in the period of deltatY=1.1+rand(1,X)*(1.3-1.1);%generate the intensity of each jumpM=prod(Y);%compute the cumulative intensity of and jumps in the period of deltat deltaw = random('Normal',0,sqrt(T),1,1); ST(s) =S0*exp((r-sigma^2/2)*T+sigma*deltaw)*M; % solution by Equation endChapter8 美式期权的定价Step 1: simulate n independent paths of the dynamics of the underlying assets Step 2: at terminal nodes, set Step 3: apply backward induction for i =T -1, …, 1 given estimated values , use regression tocalculate the coefficientSetStep 4: set 程序%Price an American put option function AmericanPut()clear allT=3; % experiation daten=3; %total number of periods S0=1;K=1.1; %exercise pricer=0.06;%risk free interest rate sigma=0.3;%vloatility of deltat=T/n;lambda=0; %Possion process with lambda per yearw = 10000; %# of replications State=[zeros(w,n),ones(w,1)]; %St=[1.00,1.09,1.08,1.34;1.00,1.16,1.26,1.54,;1.00,1.22,1.07,1.03;1.00,0.93,0.97,0.92;1.00,1.11,1.56,1.52;1.00,0.76,0.77,0.90;1.00,0.92,0.84,1.01;1.00,0.88,1.22,1.34];for s = 1:wSt(s,:) = zeros(1,n+1); St(s,1) = S0; for k = 1:nX = random('Poisson',deltat*lambda);%compute the number of jumps in the period of deltatY=1.1+rand(1,X)*(1.3-1.1);%generate the intensity of each jump M=prod(Y); %compute the cumulative intensity of and jumps in the period of deltat deltaw = random('Normal',0,sqrt(deltat),1,1); St(s,k+1) =St(s,k)*exp((r-sigma^2/2)*deltat+sigma*deltaw); % GBM endVt(s,:) = zeros(1,n+1);V(s,n+1)=max(K-St(s,n+1),0); endfor k = n:-1:2 xdata=[]; ydata=[]; for s = 1:wif St(s,k)<K % select observations for estimation of the regressionxdata=[xdata,St(s,k)];ydata=[ydata,V(s,k+1)*exp(-r*deltat)]; end endx0 = [0, 0 , 0]; % Starting guessx = lsqcurvefit(@Regfun,x0,xdata,ydata); % Call the regression function 线性回归for s = 1:w % select the decision if St(s,k)<KC(s,k)=x(1)+x(2)*St(s,k)+x(3)*St(s,k)^2; if K-St(s,k)>C(s,k); State(s,k)=1;21(2)1()()i t t Z i i S t S t e M μσ+-∆+∆+=,,ˆ()T j T T j V h X =1,ˆi jV +{},,,ˆˆmax (),()i j i i j i i jV h X C X =1110ˆˆˆn V V V n+=V(s,k)=max(K-St(s,k),0);endendendendfor s = 1:8for k = 1:n+1if State(s,k)==1;V0(s)=max((K-St(s,k)),0)*exp(-r*(k-1));breakendendendmeanV0=mean(V0)function F = Regfun(x,xdata)F = x(1) + x(2)*xdata + x(3)*xdata.^2;Chapter9 有限差分的方法美式看跌期权定价向前差分程序function imfdamput(Smax, dS, T, dT, K, R, SIG); % Implicit method% Smax : maximum stock price% dS : increment of stock price% T : maturity date% dT : time step% K : exercise price% R : risk free interest rateM = ceil(Smax/dS);%s增长的最大次数ds = Smax / M;N = ceil(T/dT); %T增长的最大次数dt = T / N;J = 1:M-1;a = .5*R*dt*J - .5*SIG^2*dt*J.^2;b = 1 + SIG^2*dt*J.^2 + R*dt;c = -.5*R*dt*J - .5*SIG^2*dt*J.^2;A = diag(b) + diag(a(2:M-1), -1) + diag(c(1:M-2),1);put = zeros(N+1, M+1);put(N+1, :) = max(K - [0:ds:Smax], 0);%赋初值put(:, 1) = K; %赋初值put(:, M+1) = 0;%赋初值for i = N:-1:1y = put(i+1, 2:M)';y(1) = y(1) - a(1)*K;put(i, 2:M) = [A \ y]';put(i, :) = max(K - [0:ds:Smax], put(i,:)); endx=0:dT:T;y=0:dS:Smax;mesh(y,x,put)%画三维图像向后差分程序f unction exfdamput(Smax, dS, T, dT, K, R, SIG); % Explicit Method% Smax : maximum stock price% dS : increment of stock price% T : maturity date% dT : time step% K : exercise price% R : risk free interest rateM = ceil(Smax/dS);ds = Smax / M; N = ceil(T/dT);dt = T / N;J = 1:M-1;a = (-.5*R*dt*J + .5*SIG^2*dt*J.^2) / (1+R*dt);b = (1 - SIG^2*dt*J.^2) / (1+R*dt);c = (.5*R*dt*J + .5*SIG^2*dt*J.^2) / (1 + R*dt);A = diag(b) + diag(a(2:M-1), -1) + diag(c(1:M-2),1);put = zeros(N+1, M+1);put(N+1, :) = max(K - [0:ds:Smax], 0);put(:, 1) = K;put(:, M+1) = 0;for i = N:-1:1y = zeros(1, M-1);y(1) = a(1)*put(i+1, 1);y(M-1) = c(M-1)*put(i+1,M+1);put(i, 2:M) = put(i+1, 2:M) * A' + y;put(i, :) = max(K - [0:ds:Smax], put(i,:)); endx=0:dT:T;y=0:dS:Smax;mesh(y,x,put)。
