主要内容

用Feynman-Kac公式模拟随机过程

本例得到了描述资产预期最终价格的偏微分方程,该资产的价格是由随机微分方程给出的随机过程。

所涉及的步骤是:

  1. 用随机微分方程定义模型参数

  2. 应用伊藤法则

  3. 解费曼-卡茨方程

  4. 计算出售资产的预期时间

1.用随机微分方程定义模型参数

一种资产价格的模型X (t)在时间间隔中定义[0, T]随机过程是由这种形式的随机微分方程定义的吗 d X μ t X d t + σ t X d B t ,在那里B (t)为单位方差参数的维纳过程。

创建符号变量T以及以下符号功能:

  • r (t)是表示即期利率的连续函数。这个利率决定了当时最终支付的折现因子T

  • u (t, x)未来价格折后的期望值是否计算为 X T 经验值 - t T r t d t 在这种情况下X(t) = X

  • μ t X 而且 σ t X 漂移和扩散是随机过程吗X (t)

信谊μ(t, X)σ(t, X)r (t) X)u (t, X)T

根据费曼-卡茨定理,u满足偏微分方程 u t + σ 2 2 2 u x 2 + μ u x - u r 0 在时间上有一个最终条件T

情商= diff (u, t) +σ^ 2 * diff (u, X, X) / 2 +μ* diff (u, X) - u * r;

数值求解器pdepe适用于初始条件。若要将最终条件转换为初始条件,请通过设置应用时间反转v(t, X) = u(t - t, X)

信谊v (t) X)eq2 = subs(eq, {u, t}, {v(t - t, X), t - t});Eq2 = feval(symengine,“改写”eq2,“差异”
eq2 =

σ T - t X 2 2 X 2 v t X 2 + μ T - t X X v t X - t v t X - v t X r T - t X

的解算器pdepe需要如下形式的偏微分方程。这里是系数cf,年代的函数xtv, v / x

c v t f x + 年代

才能解出这个方程eq2pdepe解算器,地图eq2到这个形式。首先,提取的系数和项eq2

信谊深静脉血栓形成dvXXeq3 =潜艇(eq2 {diff (v, t) diff (v, X, X)}, {dvt, dvXX});[k,terms] = coeffs(eq3, [dvt, dvXX])
k =

- 1 σ T - t X 2 2 μ T - t X X v t X - v t X r T - t X

条款=
                  
                   
                    
                     
                     
                      
                       
                        
                         
                          深静脉血栓形成
                        
                       
                       
                        
                         
                          dvXX
                        
                       
                       
                        
                         
                          1
                        
                       
                      
                     
                     
                    
                   
                  

现在,你可以重写eq2作为K (1)*terms(1) + K (2)*terms(2) + K (3)*terms(3) = 0.把带时间导数的项移到方程左边,就得到

v t k 2 2 v X 2 + k 3.

按零件手动整合右侧。结果是

X k 2 v X + k 3. - v X k 2 X

因此,要写eq2在形式上适合pdepe,使用以下参数:

C = 1;f = k(2) * diff(v, X);s = k(3) - diff(v, X) * diff(k(2), X);

2.应用伊藤法则

资产价格遵循一个乘法过程。也就是说,价格的对数可以用SDE来描述,但价格本身的期望值是有趣的,因为它描述了利润,因此我们需要一个SDE来描述后者。

一般来说,如果是随机过程X是以SDE的形式给出的,那么伊藤法则说转换后的过程G (t, X)满足

d G μ d G d X + σ 2 2 d 2 G d X 2 + d G d t d t + d G d X σ d B t

我们假设价格的对数由一维加性布朗运动给出,即:μ而且σ是常数。定义一个应用伊藤法则的函数,并用它将加性布朗运动转换为几何布朗运动。

ito = @(mu, sigma, G, t, X)...deal(mu * diff(G, X) + sigma^2/2 * diff(G, X, X) + diff(G, t),...diff(G, X) * sigma);信谊mu0sigma0[mu1, sigma1] = ito(mu0, sigma0, exp(X), t, X)
mu1 =

e X σ 0 2 2 + μ 0 e X

sigma1 =
                  
                   
                    
                     
                      
                       
                        
                         σ
                       
                       
                        
                         0
                       
                      
                      
                      
                      
                       
                        
                         e
                       
                       
                        
                         X
                       
                      
                     
                    
                   
                  

取代exp (X)通过Y

信谊Ymu1 = subs(mu1, X, log(Y));sigma1 = subs(sigma1, X, log(Y));f = f(t, log(Y));s = s(t, log(Y));

为简单起见,假设利率为零。这是一种特殊情况也称为Kolmogorov反方程。

R0 = 0;C = subs(C, {mu, sigma}, {mu1, sigma1});S = subs(S, {mu, sigma, r}, {mu1, sigma1, r0});F = subs(F, {mu, sigma}, {mu1, sigma1});

3.解费曼-卡茨方程

在将符号表达式转换为MATLAB函数句柄之前,必须替换函数调用,例如(v(t, X), X)而且v (t) X),有变量。您可以使用任何有效的MATLAB变量名。

信谊dvdxV;dvX = diff(v, X);c = subs(c, {dvX, v}, {dvdx, v});f = subs(f, {dvX, v}, {dvdx, v});s = subs(s, {dvX, v}, {dvdx, v});

对于平面几何(平移对称),设置以下值:

M = 0;

另外,将数值赋给符号参数。

Muvalue = 0;σ值= 1;C0 = subs(c, {mu0, sigma0}, {muvalue, sigmavalue});F0 = subs(f, {mu0, sigma0}, {muvalue, sigmavalue});S0 = subs(s, {mu0, sigma0}, {muvalue, sigmavalue});

使用matlabFunction创建一个函数句柄。传递系数c0f0,s0所要求的格式pdepe,即带有三个输出参数的函数句柄。

FeynmanKacPde = matlabFunction(c0, f0, s0,“var”, [t, Y, V, dvdx])
FeynmanKacPde =Function_handle with value:@ (t Y V, dvdx)协议(1.0,(Y ^ 2 * dvdx)。/ 2.0,(Y * dvdx) / 2.0)

作为最后的条件,取恒等映射。也就是时刻的收益 T 由资产价格本身给出。你可以修改这条线来研究衍生工具。

费曼卡奇= @(x) x;

偏微分方程的数值求解只能应用于有限域。因此,必须指定边界条件。假设资产在其价格高于或低于某一限制的时刻被出售,因此解决方案v必须满足X - v = 0在边界点x.你可以选择另一个边界条件,例如,你可以使用V = 0模拟淘汰选项。第二个和第四个输出中的零表示边界条件不依赖于 d v d x

FeynmanKacBC = @(xl, vl, xr, vr, t)...Deal (xl - vl, 0, xr - vr, 0);

选择空间网格,这是价格的值范围x.将左侧边界设为零:它不能通过几何随机游走到达。将正确的边界设置为1:当资产价格达到此限制时,资产将被出售。

Xgrid = linspace(0,1,101);

选择时间网格。因为在开始时应用了时间反转,所以它表示到最后时间的剩余时间T

Tgrid = linspace(0,1,21);

解这个方程。

sol = pdepe(m,FeynmanKacPde,FeynmanKacIC,FeynmanKacBC,xgrid,tgrid);

画出解。预期销售价格几乎与当时的价格成线性关系t,也弱上t

Surf (xgrid, tgrid, sol)“资产的预期售价”);包含(“价格”);ylabel (“时间”);zlabel (“预期最终价格”);

图中包含一个轴对象。标题为“资产预期售价”的坐标轴对象包含一个类型为surface的对象。

带漂移的几何布朗运动的状态 μ 1 在时间t对数正态分布随机变量是否具有期望值 经验值 μ 1 t 乘以它的初值。它描述了一种资产的预期售价,这种资产因为达到了一个极限而永远不会被出售。

预期=转置(exp(tgrid./2)) * xgrid;

用上面得到的解除以这个期望值,可以看出在限价时过早卖出会损失多少利润。

sol2 = sol./期望;Surf (xgrid, tgrid, sol2)'预期最终价格之比:有/没有销售订单x=1')包含(“价格”);ylabel (“时间”);zlabel (“最终价格比率”);

图中包含一个轴对象。带有标题的坐标轴对象包含一个类型为surface的对象:期望最终价格的比率:在x=1处有/没有销售订单。

例如,绘制在一段时间内已下达限价销售订单的资产与未下达限价销售订单的同一资产的预期收益之比T的函数t.考虑订单限额分别为当前价格的两倍和四倍的情况。

Plot (tgrid, sol2(:, xgrid == 1/2));持有Plot (tgrid, sol2(:, xgrid == 1/4));传奇(“限价为现价的两倍”“限价为现价的四倍”);包含(“订单有效的时间跨度”);ylabel (“最终价格比率:有/没有销售订单”);持有

图中包含一个轴对象。axis对象包含2个line类型的对象。这些标的代表两倍现价的限价,四倍现价的限价。

4.计算出售资产的预期时间

这是一个教科书式的结果,当达到限制并出售资产时,预期退出时间由下式给出:

信谊y (X)exitTimeEquation (X) =潜艇(eq, {r u}, {0 y (X)}) = = 1
exitTimeEquation (X) =

σ t X 2 2 X 2 y X 2 + μ t X X y X - 1

此外,y边界处必须为零。为μ而且σ,插入实际考虑的随机过程:

exitTimeGBM = subs(subs(exitTimeEquation, {mu, sigma}, {mu1, sigma1}), Y, X)
exitTimeGBM (X) =

X σ 0 2 2 + X μ 0 X y X + X 2 σ 0 2 2 X 2 y X 2 - 1

求解任意区间边界的问题一个而且b

信谊一个bexitTimeGBM, y(a) == 0, y(b) == 0)
exitTime =

2 μ 0 σ 5 日志 b - 2 μ 0 σ 4 日志 一个 + 一个 σ 1 σ 0 2 σ 5 σ 4 - b σ 1 σ 0 2 σ 5 σ 4 σ 3. - 日志 X μ 0 + σ 2 σ 7 - σ 6 - 一个 σ 1 σ 0 2 σ 5 + b σ 1 σ 0 2 σ 4 σ 3. + X σ 1 σ 0 2 σ 2 2 μ 0 2 在哪里 σ 1 2 μ 0 σ 0 2 σ 2 e - 2 μ 0 日志 X σ 0 2 σ 3. 2 μ 0 2 σ 5 - σ 4 σ 4 e - σ 6 σ 0 2 σ 5 e - σ 7 σ 0 2 σ 6 2 μ 0 日志 b σ 7 2 μ 0 日志 一个

因为你不能代替Mu0 = 0直接计算极限Mu0 -> 0代替。

L = limit(subs(exitTime, sigma0, sigmavalue), mu0, muvalue)
L =
                  
                   
                    
                     
                      
                       -
                      
                       
                        
                         
                          
                           
                            
                             
                              
                               日志
                             
                             
                              
                              
                               
                                
                                 X
                               
                              
                              
                             
                            
                            
                             -
                            
                             
                              
                               日志
                             
                             
                              
                              
                               
                                
                                 一个
                               
                              
                              
                             
                            
                           
                          
                         
                        
                        
                        
                        
                         
                          
                           
                            
                             
                              
                               日志
                             
                             
                              
                              
                               
                                
                                 X
                               
                              
                              
                             
                            
                            
                             -
                            
                             
                              
                               日志
                             
                             
                              
                              
                               
                                
                                 b
                               
                              
                              
                             
                            
                           
                          
                         
                        
                       
                      
                     
                    
                   
                  

l未定义A = 0.假设0 < x < 1

假设(0 < X & X < 1)

使用值B = 1对于右边的边界,计算极限。

极限(下标(L, b, 1), a, 0,“对”
ans =
                  
                   
                    
                   
                  

预计退出时间是无限的!