随机存储问题的模拟方法：一种统计方法

（10） 从概念上讲，我们有一个从值函数V到所选m的映射*和c*, 编码asm*: （t，P，I，m）7→ J、 连续值q、值函数V和控制映射m的显式依赖性*在制度上，m是转换成本K（mtk，mtk+1）的结果，其缺席将导致需要解决的低维问题，见备注3。相反，方程式（6）提供了价值函数的表示，作为基于最优政策m的未来收益的条件预期*. 因此，任何估计值^m：（t，P，I，m）7→ Jof控制图自然地导出值函数的相应估计值^V。具体而言，^m产生动力学^Itk+1=^Itk+actk（^mtk+1（t，P，I，^mtk））t、 进而可用于样本外正向模拟，^V（0，P，I，m）=E“K-1Xs=1e-rtsπ（Pts，^mts，^mts+1）+e-rTW（PT，^IT）#。虽然上面没有明确显示^Its，但它是推动^mts+1的关键因素（t、Pts、^Its、^mts）。图1显示了天然气储存设施的（Pt）轨迹和（^It）的几个对应部门，并通过其初始库存I（视为外部参数）对其进行了索引，从而说明了这种双重联系（更多详细信息，请参见第6.1节）。一个有趣的观察结果是，时间t库存对I=I的依赖性相当弱，即库存水平在初始“瞬态”时间段后合并：时间t库存=I。请注意，由于控件是指定的反馈形式，一旦^Iit=^Iitwe有^mits=^mits，库存路径将永远保持在一起。该图还说明了“低买高卖”的基本准则：当Ptis低时，控制库存^为高（并不断增加），当Ptis高时，^为低（并不断缩小）。因此，我们看到围绕最小和最大存储级别Imin、Imax的^Ita聚类，这表明有限存储容量施加的强大约束。

为了可视化估计的最佳价格0.5 1 1.5 2 2.5 32468100 0.5 1.5 2 2.5 3时间（年）050100015002000库存图1：顶部面板：在（24）中对数均值回归过程后商品价格（Pt）的给定轨迹。下面板：受控库存的相应轨迹^Itstartingat^I∈ {0, 500, 1000, 2000}. 该图与第6.1节使用PR-1D解决方案的储气库示例相关。政策^m，图2a绘制了控制图（P，I）7→ ^m（t，P，I）在固定的时间步长t，即- 对于第6.1节中的储气示例，t=0.3年。状态空间分为三个区域：当P较高时，最佳提取：^m=-1.如果P较低，注入^m=+1为最佳；在中间，或者如果库存非常大，最好什么都不做^m=0。通常，该图是通过计算当前库存I并将^m作为P的函数来解释的。然后，我们可以根据注入/退出边界Binj（I，t）和Bwdr（I，t）：Binj（I，t）：=sup{Pt:^m（t，Pt，I）=+1}，Bwdr（I，t）：=inf{Pt:^m（t，Pt，I）=-1}. （11） 由于注入价格较低，Binj（I，t）代表最高价格，而注入（^m（t，Pt，I）=+1）是最佳策略。类似地，Bwdr（I，t）表示最低价格，其中提取（^m（t，Pt，I）=-1） 是最佳策略。区间[Binj（I，t），Bwdr（I，t）]是无作用区。这些边界绘制为T的函数- 图2b中三种不同库存水平下的t。

一个突出的特征是边界层- t型→ 0，其中政策主要由终端罚款W（P，I）而非直接利润因素驱动。在本例中，由于“曲棍球棒”惩罚W（PT，IT），迫使控制器以库存水平I=1000为目标，如T- t型→ 0时，注射成为I<1000的最佳策略，与价格无关，而退出则成为I>1000的最佳策略（无行动区域有效消失）。相反，对于大T-t、 边界t 7→ Binj（I，t）和t 7→ Bwdr（I，t）本质上是时间平稳的。最后，图2c显示了价值函数^V（0，P，I）作为价格和库存坐标中的二维曲面。正如之前的研究所指出的，对于固定价格P，我们观察到价值和库存之间存在线性关系。然而，作为P的函数，V（0，·，I）在低库存水平下呈非线性递减，在大库存水平下呈非线性递增。（a） 对照map0 0.05 0.3 0.45 0.6 0.75 0.9 1.05 1.2到期时间（年）23456789101112定价WDR（1000，t）Binj（1000，t）Bwdr（500，t）Binj（500，t）Bwdr（1500，t）Binj（1500，t）（b）产品/注入。边界（c）值函数图2：左面板：t=2.7年时控制图^m（t，P，I）的快照。中间：注射BINJ（I，·）和退出Bwdr（I，·）边界作为I的时间函数∈ {500, 1000, 1500}.右：t=0时的值函数^V（0，P，I）。结果通过常规设计和PR-1D回归得到。3动态仿真算法我们考虑的数值算法为（10）中的连续值提供近似值，表示为^q（·）。然后，我们通过反向归纳ontk得到V（·）的递归解。即，从已知的终端条件v（T，PT，IT，m）=q（tK，PtK，ItK，m）=W（PT，IT）开始m、 我们应用反向归纳法估计^q（tk，·，·，·，·），作为k=k- 1到k=0。

递归构造确保在步骤tkwe知道连续函数^q（tk+1，·，·，·，·），因此可以找到^m（tk+1，P，I，m），如（9）所示。然后取^V（tk，P，I，m）=π（P，m，^mtk+1）+^q（tk+1，P，I+a（c（^mtk））t、 ^mtk+1）（12）并通过学习^q（tk，·）作为e的条件期望来完成归纳步骤-rt^V（tk+1，·）。备注2。在上面，我们将连续值视为一步aheadvalue函数的条件期望，与原始的[35]方案相匹配。更一般地，我们可以使用最佳制度选择m展开动态规划方程（8）*tk+1和相应的控制装置c*TK适用于任何w≥ 1 asV（tk，Ptk，Itk，mtk）=E[π（Ptk、mtk、m*tk+1）+q（tk，Ptk，Itk+1，m*tk+1）Ptk]=Ehπ（Ptk、mtk、m*tk+1）+e-rtπ（Ptk+1，m*tk+1，m*tk+2）+e-rtq（tk+1，Ptk+1，I*tk+2，m*tk+2）Ptki…=Ehvk：k+w（π，q）（Ptk，Itk，mtk）Ptk、Itk、m*tki（13）在路径增益Vk方面：k+w（π，q）（Ptk，Itk，mtk）：=k+w-1Xs=ke-r（s）-k）tπ（Pts、mts、m*ts+1）+e-rw公司tq（tk+w、Ptk+w、Itk+w+1、mtk+w+1）。（14） 类似地，连续值q（tk，Ptk，Itk+1，mtk+1）可以写为q（tk，Ptk，Itk+1，mtk+1）=Ehvk+1：k+w（π，q）（Ptk+1，Itk+1，mtk+1）| Ptk，Itk，mtki。（15） （15）中的部分路径构造可以追溯到[15，16]。它包括西西里斯·范罗伊（TvR）算法（其中w=1）和朗斯塔夫·施瓦茨（CLS）算法（其中w=K）- k- 1、在本文的其余部分，我们将x w=1，并使用vk的短符号：k+1（π，q）（Ptk，Itk，mtk）作为vk（Ptk，Itk，mtk）或仅作为vk。基于模拟的框架依赖于生成路径特性vk+1（Ptk+1，Itk+1，mtk+1），并通过蒙特卡罗近似将连续函数^q（tk，·）恢复为（15）中的预期值。也就是说，我们的目标是使用回归过程将vk+1投影到近似空间HK上：ˇq（tk，·）：=arg minhtk∈Hkkhtk公司- vk+1k。

（16） 作为规范设置，Hk=span（φ，…，φR）是由基函数φi生成的线性空间，因此近似值ˇq（tk，·）=PRi=1βiφi（·）通过其系数向量β来描述。为了估计β，我们基于模拟设计Dk：=（Pntk，Intk+1）Nn=1和相应的实现路径值vnk+1，n=1，…，来解决离散优化问题，N从（Pntk，Intk+1，mtk+1）开始的轨迹：^q（tk，·，·，mtk+1）=arg minhtk∈HkNXn=1 | htk（Pntk，Intk+1）-vntk+1 |。（17） 因此，^q是投影^q的经验近似值，具有相应的有限样本误差。特别是，根据Dk上的温和条件，我们有^q→ _q.传统上（17）是通过为外部（价格）过程生成N条全局路径（Pnt）来实现的，该路径从t=0开始直到成熟度t，然后永久存储在内存中，用于整个反向归纳，引入了大量开销。算法1将其替换为设计DK和相关的单步轨迹（Pnt:t+1，Int+1）。备注3。在没有转换成本K（i，j）的情况下≡ 0i、 回归问题的维数可以从3降到2。事实上，由于库存过程受到完全控制，持续值仅取决于库存和制度（Itk，mtk），通过下一步库存Itk+1=Itk+a（c（mtk））t、 类似地，值函数与当前状态V（tk、Ptk、Itk）无关。然后，我们使用（Ptk，Itk+1）生成的投影子空间编写EQ（tk，Ptk，Itk+1）。当存在转换成本时，在回归过程中可能会出现相同的减少，但当前制度MtKre是状态的一部分，因为它会影响延续值q，因此不会实现总体维度节约。

实际上，这是通过为每个∈ J如（17）所示。为了应用（17），我们还需要计算任意未来状态（PNT、Ints、MNT）下的连续值，s>k对应于预测之前的^q（ts，·，·，·，·）。这两种操作拟合和预测是统计近似程序的主要工作。上述程序生成估计的连续值^q（tk，·，·，·，·，·），对应的最佳状态^m（tk+1，·，·，·，·，·）匹配（9），从而生成控制c（^m）。最后，通过生成样本路径（Pn0:t，^In0:t，mn0:t），t=0，P，I，m时的值函数V（0，P，I，m）在样本外不近似，其中路径库存^I基于刚刚估计的最优控制图^m，参见图1。因此，^V（0，P，I，m）=NPNn=1v0:T（π，^q）（0，Pn，In，^mn），其中v0:T（π，^q）（0，Pn，In，^mn）是使用每个样本路径的^m估计的总累计贴现利润。由于策略^m必然是次优的，因此^V（0，·）是真V的下界，模化最后一次平均中使用的nTrajelector的蒙特卡罗误差。备注4。请注意，最终估计值在形式上是样本外模拟Pn0:T和样本内模拟（Pnt:T+1）的函数。为了便于比较不同的方法，我们尽可能使用“测试”情景数据库（Pn0:T），通过该数据库，我们可以直接评估在价格过程的给定样本路径上获得的不同控制/累计收入。动态仿真算法1（DEA）提供了解决存储问题的总体模板。

之所以选择该术语，是因为该算法允许用户随时间改变投影空间HK和模拟设计DK，因此它是“动态的”；以及“模拟”，以反映统计学习视角，其目标是近似或模拟连续值q（t，·）。其总体复杂性为O（KN）。相对于现有文献，DEA增加了以下工具：o可能包括跨时间步tk的不同回归空间HK；o一般模拟设计Dk，同样可能随时间步而变化；o消除了在内存中存储全局价格路径的要求，而是“在线”模拟新路径，即在每个时间步。这还允许在时间步长tk改变模拟数nk。（17）的另一种选择是在投影期间用作状态变量（Ptk+1，Itk+1），然后算法1：动态仿真算法（DEA）-O（KN）数据：K（时间步），（Nk）（每步仿真预算）1生成设计DK-1，m：=（PDK-1，mK-1，IDK-1，mK），尺寸为NK-1每米∈ J、 2生成一步路径Pn、DK-1，mK-17→ Pn，DK-1，对于n=1，…，mk，NK公司-1和m∈ J3终端条件：vnK，m← W（Pn，DK-1，mK，In，DK-1，mK）对于n=1，NK公司-1和m∈ k=k时为J4- 1.m为1 do5∈ J do6^q（k，·，·，m）← arg minhk公司∈HkPNkn=1 | hk（Pn、Dk、mk、In、Dk、mk+1）-vnk+1，m | 7生成设计Dk-1，m：=（PDk-1，mk-1，IDk-1，mk），尺寸为Nk-1每米∈ J8生成一步路径Pn、Dk-1，mk-17→ Pn，Dk-1，对于n=1，…，mk，Nk公司-19 end10表示n=1。

，Nk-1和m∈ J do11米← arg最大值∈J{π（Pn，Dk-1，mk，m，j）+^q（k，Pn，Dk-1，mk，In，Dk-1，mk+a（ck（j））t、 j）}12 vnk，m← π（Pn，Dk-1，mk，m，m）+e-rt^q（k，Pn，Dk-1，mk，In，Dk-1，mk+a（ck（m））t、 m）13端14端15返回{q（k，·，·，m）}k-1k=1，m∈j分析条件期望：^q（tk，P，I，m）=E“arg minhtk+1∈Hk+1NXn=1 | htk+1（Pntk+1，Intk+1）-vnk+1|Ptk=P，Itk+1=I#。（18） 这被称为后回归蒙特卡罗（RLMC），并降低了估计的^q中的方差【31】。然而，对条件期望的闭式表达式的要求通常不使用RLMC和非参数近似空间。以下两个部分提供了DEA中两个主要步骤的菜单：选择近似空间HK和设计Dk。注意，这些可以在tk中混合和匹配。总的来说，这两个步骤对应于一个机器学习任务：给定一个随机模拟器（在给定初始条件（t，P，I，m）的情况下返回一步前的结果），我们希望学习输入输出关系，即预测任何（样本内或新的样本外）输入的预期响应。因此，在该语言中，DEA是一个（递归的）学习任务序列，性能度量由最终答案的质量（0，·，·，·，·，·）给出。最后，我们在本节结束时对算法的另一个属性进行了评论，即前瞻参数wwhich，我们在这项工作中没有探讨它。4近似空间在本节中，我们确定给定的时间步长tk，并考虑将连续值q（tk，P，I，m）近似为函数htk（P，I）的问题。下面我们一般使用x=（Pntk，Intk+1）Nn=1和y=（vnk+1）Nn=1来表示回归过程中使用的数据集。我们注意到，虽然完整的状态空间通常是（P，I，m），但由于m是一个因子变量（而不是r值），它是离散处理的，即。

对于m的每个级别，构造一个单独的近似值^h（·，m）∈ J（因此，在典型存储设置中有三个单独的近似值）。当不存在转换成本且连续值的维数仅为2时，后一个问题就没有意义了。统计假设是，投入产出关系由yn=h（xn）+σξ描述，其中ξ~ N（0，1），（19），其中h∈ hk是要学习的未知函数，σξ是噪声。在我们的案例中，噪声是由于（Ptk+1）相对于Ptk的随机冲击引起的，这导致实现的路径效应相对于平均连续值的变化。经典回归框架是一个线性模型，其中hk是由一些基函数（φi）跨越的向量空间。然后，在一般x*由回归系数β：h（x）控制*) =~βT~φ（x*). 然后，关键的挑战是指定基函数φi，因为固有近似误差（即真实q（tk，·）和Hk之间的距离）强烈影响解的质量。因此，文献中探索了各种回归方法，例如，全局多项式回归[4，7，10]，径向基函数[30]，支持向量回归[29]，核回归[26]，神经网络[19]和分段线性回归[37]。受最近百慕大期权研究[28]的启发，下面我们还介绍了高斯过程（GP）回归在解决存储问题中的应用。据我们所知，我们是第一篇在这种情况下使用GPs的论文。4.1二元分段近似继Longstaff-Schwartz关于百慕大期权定价的开创性工作【27】之后，研究人员还将多项式回归（PR）用于存储问题【4、7、10】。

r次全局多项式逼近htk=Piβiφi，具有（r+1）（r/2+1）基函数，取φi（P，i）=Pα（i）iα（i），其中基函数的总阶数为α+α≤ r、 例如，全局二次近似（r=2）有6个基函数{1，P，P，I，I，P·I}，而三次PR有10个基函数。我们的实验表明，由于对延拓值形状的严格限制以及由此产生的误差反向传播，PR通常会导致较差的性能。一种流行的替代方法是使用基于（P，I）空间划分的分段近似，限制为[min1≤n≤NPnti，最大值1≤n≤NPnti]×[Imin，Imax]，进入M=MP×多角子域Di，i，i=1，2，MP；i=1，2，惯性矩。然后，我们考虑{φi，ig}形式的基函数，支持度限制为Di，i。例如，我们可以采用分段线性近似，即g=1，2，3，φi，i（P，i）=1（P，i）∈Di，iφi，i（P，i）=P·1（P，i）∈Di，iφi，i（P，i）=i·1（P，i）∈总的来说，我们需要估计3MP系数。还可以添加更高阶的项，例如交叉项P·I或二次项P，I。分段回归提供了一种分而治之的优势，通过跨Di，I的循环进行整体拟合；在每种情况下，仅选择一小部分数据来学习一些系数。这减少了回归子步的总体工作量，并允许并行处理。相对于PR，分段回归对拟合连续值的任意形状也更“稳健”。它们的主要缺点是子域边界上^q的内在不连续性，需要指定MP，mian，然后构造矩形子域Di，i。