金融工程与并行计算：第三章 多资产模拟与R的使用 Part 2

第三节 R的Cholesky函数

在多资产的模拟中，我们需要将相关性矩阵进行Cholesky分解。形成下三角与上三角矩阵的乘积。以一个2╳2的相相关性矩阵为例。M可分解成L与U的乘积。

...........................................................(3.3.1)

 

#001 R version 3.1.1 (2014-07-10) --"Sock it to Me"

#002 Copyright (C) 2014 The R Foundation for Statistical Computing

#003 用 'demo()' 来看一些示范程序，用 'help()' 来检视线上辅助档案，或

#004 用 'help.start()' 透过 HTML 浏览器来看辅助档案。

#005 用 'q()' 离开 R。

#006

#007 > M <- matrix(c(1.0, 0.8, 0.8,1.0), nrow = 2, ncol = 2, byrow = TRUE)

#008 > U <- chol(M)

#009 > U

#010     [,1] [,2]

#011 [1,]   1  0.8

#012 [2,]   0  0.6

#013 > L <- t(U)

#014 > L

#015     [,1] [,2]

#016 [1,] 1.0  0.0

#017 [2,] 0.8  0.6

#018 > A = L %*% U

#019 > A

#020     [,1] [,2]

#021 [1,] 1.0  0.8

#022 [2,] 0.8  1.0

#023 >

程序行表3.1

 

程序行表3.1示范了R语言的Cholesky分解函数，#007使用matrix()函数建立M矩阵；先以c()函数建立一维向量，内含四个分量，当作传入数据；然后，说明M矩阵为两列，两行的矩阵，并以列优先的方式，读取c()函数建立的一维向量。#008使用chol()函数执行Cholesky分解，将上三角矩阵传出给U矩阵。#013使用t()函数进行U矩阵的转置，将下三角矩阵传出给L矩阵。#018以L %*% U进行矩阵相乘，结果传给A。如果想知道每个变量的内容，只要在提示符号后输入变量的名称，按Enter就可看到结果。如#009输入U按Enter，便可看到2╳2的上三角矩阵。

第四节 C#多资产的模拟

假设有两资产，契约为买权型式，偿付为两者的价差减去Strike。

...........................................................................................(3.4.1)

契约为一年期，模拟的架构同前。

#001 namespace STBiBSProcess

#002 {

#003     public partial class Form1 : Form

#004     {

#005         publicForm1()

#006         {

#007             InitializeComponent();

#008         }

#009         privatevoid button1_Click(objectsender, EventArgs e)

#010         {

#011             Application.Exit();

#012         }

#013        

#014         const int BlocksPerGrid = 256;

#015         const int ThreadsPerBlock = 256;

#016         static int NPath;

#017         static int MStep;

#018         static DateTime RefDate;

#019         static List<DateTime> FixingDate;

#020         static int[] h_StepGrid;

#021

#022         privatevoid button2_Click(objectsender, EventArgs e)

#023         {

#024             doubleAsset1 = 100.0;

#025             doubleAsset2 = 100.0;

#026             doubleStrike = 1.0;

#027             double TTM =1.0;

#028             doubleSigma = 0.3;

#029             doubleRate = 0.04;

#030             doubleYield = 0.02;

#031             doubledt = 1.0 / 365.0;

#032             doubleRho = 0.5;

#033             doubleSqrt1mRho2 = Math.Sqrt(1 - Rho * Rho);

#034

#035             //RandomRnd = new Random(1234);

#036             MersenneTwisterRnd = new MersenneTwister(1234);

#037             StopwatchSW = new Stopwatch();

#038             NPath = ThreadsPerBlock *BlocksPerGrid;

#039             RefDate = newDateTime(2014, 7, 1);

#040             FixingDate = new List<DateTime>();

#041

#042             //MStepdepend on the remaining FixingDay Number, Including RefDate

#043             MStep = 13; // 0, 1, 2,..., 12, 0 for RefDate

#044             for(int i = 0; i < MStep; i++)

#045             {

#046                 FixingDate.Add(RefDate.AddMonths(i));

#047             }

#048

#049             h_StepGrid = new int[MStep];

#050             for(int i = 0; i < MStep; i++)

#051             {

#052                 TimeSpanTS = FixingDate[i].Subtract(RefDate);

#053                 h_StepGrid[i] = (int)TS.TotalDays;

#054             }

#055

#056             double[,]S1 = new double[NPath,MStep];

#057             double[,]S2 = new double[NPath,MStep];

#058             double[]Value = new double[NPath];

#059             doubles1, s2, n1, n2, e1, e2;

#060             SW.Start();

#061             for(int i = 0; i < NPath; i++)

#062             {

#063                 s1 = Asset1;

#064                 s2 = Asset2;

#065                 intdiff = 0;

#066                 S1[i, 0] = s1;

#067                 S2[i, 0] = s2;

#068                 for(int j = 0; j < (MStep - 1); j++)

#069                 {

#070                     diff = h_StepGrid[j + 1] -h_StepGrid[j];

#071                     for(int k = 0; k < diff; k++)

#072                     {

#073                         e1 = DStat.N_Inv(Rnd.NextDouble());

#074                         e2 = DStat.N_Inv(Rnd.NextDouble());

#075                         n1 = e1;

#076                         n2 = Rho * e1 +Sqrt1mRho2 * e2;

#077                         s1 = s1 * Math.Exp(((Rate-Yield) - (Sigma*Sigma)/2.0) * dt

#078                             + (Sigma * Math.Sqrt(dt) * n1));

#079                         s2 = s2 * Math.Exp(((Rate-Yield)- (Sigma*Sigma)/2.0) * dt

#080                             + (Sigma * Math.Sqrt(dt) * n2));

#081                     }

#082                     S1[i, j + 1] = s1;

#083                     S2[i, j + 1] = s2;

#084                 }

#085                 Value[i] = Math.Max((s1-s2) - Strike, 0);

#086             }

#087

#088             doublesum = 0.0;

#089             for(int i = 0; i < NPath; i++)

#090             {

#091                 sum = sum + Value[i];

#092             }

#093             sum = (sum / NPath) * Math.Exp(-(Rate - Yield) * TTM);

#094             SW.Stop();

#095                    

#096             textBox1.Text =SW.ElapsedMilliseconds.ToString();

#097             textBox2.Text = sum.ToString();

#098             for(int j = 0; j < MStep; j++)

#099             {

#100                 listBox1.Items.Add(S1[0,j].ToString());

#101                 listBox2.Items.Add(S2[0,j].ToString());

#102             }

#103         }

#104     }

#105 }

程序行表4.1

#073与#074产生两个独立的标准常态变量，#075与#076使用(3.1.2)产生两个具相关性的随机变量，#077与#078以之产生两个具相关性的价格程序。执行结果如下。

 

 

读者可以调整不同的相关性，看看此参数对价格的影响。事实上，这是一个相当敏感的参数，忽视或估计错误，对结果有很大的影响。