金融工程与并行计算：第三章 多资产模拟与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以之产生两个具相关性的价格程序。执行结果如下。
读者可以调整不同的相关性，看看此参数对价格的影响。事实上，这是一个相当敏感的参数，忽视或估计错误，对结果有很大的影响。