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

第五节 R的矩阵运算

R语言提供了很多矩阵的运算功能。

一、向量的基本操作

绝大部分R的数据(包括常数与变量)都以内建向量的形式呈现，这种内建向量称之为原型向量，维度为一。单一常数也是一种原型向量，只是其长度为一的一维向量。在本书中，一个空向量我们以[]表示，[1,2, 3]表示一个有三个整数元素的向量。

(一)初始化

A.空向量

产生一个长度为0的空值的实数向量，

VecR = numeric()

产生一个长度为0的空值的整数向量，

VecI = integer()

B.长度为n的向量

产生一个长度为n，初值为零的实数向量，足标由1开始，1…n，

VecR = numeric(n)，VR =numeric(10)

产生一个长度为n，初值为零的整数向量，

VecI = integer(n)，VI =integer(10)

C.设定初值的向量

产生一个长度为n，初值为1到n的整数向量，

VecI = 1:n，VI = 1:30

产生一个初值为n到m，递增u的整数向量，

VecI = seq(from = n, to = m, by = u)，VI= seq(from = -5, to = 5, by = 2)

产生一个长度为n，初值为x₁、x₂、…、x_n的实数向量，

VecR = c(x₁, x₂,…, x_n)，VR= c(1.1, 2.3, 3.7)

产生一个初值为n到m，递增u的实数向量，

VecR = seq(from = n, to = m, by = u)，VR= seq(from=1.0, to=2.5, by=0.1)

(二)型别判断与转换

判断资料型别，是否实数，is.numeric()。Flag= is.numeric(VecR)，Flag = TRUE。

判断资料型别，是否整数，is.integer()。Flag= is.integer(VecR)，Flag = FALSE。

资料转型为实数，as.numeric()。VecR =as.numeric(VecI)。

资料转型为整数，as.integer()。VecI =as.integer(VecR)。

(三)元素存取

A.单一元素

Z = numeric(10)；Z[i], i=1…10

Z[1] = 1.1；Z[2] = 2.3；…

B.区段元素

Z[1:3] = c(1.1, 2.3, 3.7)

X = Z[c(1,3)]，X = c(Z[1],Z[3]) = [1.1, 3.7]

C.删除元素

Z[1:3] = c(1.1, 2.3, 3.7)，Z[-1]，Z= c(Z[2], Z[3])

Z[1:3] = c(1.1, 2.3, 3.7)，Z[-c(2,3)]，Z = c(Z[1])

(四)基本运算

一般运算是针对每一元素运作，

A.函数对每一元素运作

abs(X)

sqrt(X)

log(X)

exp(X)

sin(X)

ceiling(X)

floor(X)

trunk(X)

round(X, digits = s)

 

B.函数对整个向量运作

length(X)：向量长度

min(X)

max(X)

mean(X)

sum(X)

prod(X)

sd(X)

var(X)

median(X)

C.运算规范

X = (1, 0, -1, 2)，Y = (-2, 1,0, -2)，Z = X * Y，Z = (-2, 0, 0, -4)，内积。

X = (1, 0, -1, 2)，Y = 2 * X =(2, 2, 2, 2) * X = (2, 0, -2, 4)，长度较短者，循环使用

X = (1, 0, -1, 2)，Y = X * (1:2)= X * (1, 2, 1, 2) = (1, 0, -1, 4)，长度较短者，循环使用

二、矩阵的基本操作

(一)初始化

空矩阵，m列n行。

MatA = matrix(nrow = m, ncol = n)

初值为零矩阵，m列n行。

MatA = matrix(0, nrow = m, ncol = n)

单位矩阵，m列m行。

MatI = diag(m)

初值1~12矩阵，4列3行，列优先。

MatA = matrix(1:12, nrow = 4, ncol = 3,byrow = TRUE)

初值1~12矩阵，4列3行，行优先。

MatA = matrix(1:12, nrow = 4, ncol = 3)

 

(二)元素存取

初值1~12矩阵，4列3行，列优先。

MatA = matrix(1:12, nrow = 4, ncol = 3,byrow = TRUE)

MatA =

1              2              3              4

5              6              7              8

9              10            11            12

A23 = MatA[2, 3]

A23 = MatA[8]      #bycol index

取整列，R1为向量

R1 = MatA[1,]

取整行，C1为向量

C1 = MatA[,1]

 

(三)合并运算

行结合(列相同)

MatC = cbind(MatA, MatB)

列结合(行相同)

MatC = rbind(MatA, MatB)

 

(四)基本运算

矩阵加、减，

MatA + MatB

MatA - MatB

矩阵相乘，

MatA %*% MatB

矩阵点乘，

MatA * MatB = MatA[i, j] * MatB[i, j]

矩阵转置，

TMatA = t(MatA)

反矩阵，

InvMatA = solve(MatA)

矩阵行列式值，

D = det(MatA)

 

三、矩阵操作函数

(一)方程式求解

MatA * Vx = Vb

Vx = solve(MatA, Vb)

 

(二)矩阵之光谱分解

EV = eigen(MatA)

EigenValues = EV$val

EigenVectors = EV$vec

 

(三)矩阵之QR分解

MatM = MatU %*% MatD %*% t(MatV)

SVD = svd(MatM)

 

(四)矩阵之奇异值分解

Vy = MatX * Vb

QR = qr(MatX)

b = QR.coef(MatX, Vy)

fit = QR.fitted(MatX, Vy)

res = QR.resid(MatX, Vy)

 

四、多维数组的操作

(一)宣告

以一个三维空数组为例，

> A <-array(dim= c(2,3,4))

> A

, , 1

     [,1] [,2] [,3]

[1,]   NA  NA   NA

[2,]   NA  NA   NA

 

, , 2

     [,1] [,2] [,3]

[1,]   NA  NA   NA

[2,]   NA  NA   NA

 

, , 3

     [,1] [,2] [,3]

[1,]   NA  NA   NA

[2,]   NA  NA   NA

 

, , 4

     [,1] [,2] [,3]

[1,]   NA  NA   NA

[2,]   NA  NA   NA

 

表示A为三维(x,y, z)数组，x轴有2个成分，y轴有3个成分，z轴有4个成分。

 

(二)设定初值

如果要设定初值，可以如下设定，

> B <- array(1:24,dim= c(2,3,4))

> B

, , 1

     [,1] [,2] [,3]

[1,]    1    3    5

[2,]    2    4    6

, , 2

     [,1] [,2] [,3]

[1,]    7    9   11

[2,]    8   10   12

, , 3

     [,1] [,2] [,3]

[1,]   13   15   17

[2,]   14   16   18

, , 4

     [,1] [,2] [,3]

[1,]   19   21   23

[2,]   20   22   24

 

注意初值的设定顺序，(x,y)二维矩阵是bycolumn设定，每一(x, y)二维矩阵完成后才会设定下一个二维矩阵。

 

(三)操作

数组以其足标存取其值，

> B[1, 1, 1] = 51

> Fst = B[1, 1, 1]

> Fst

[1] 51

 

 

第六节 绘图功能

R语言有很好的绘图功能，下面把第四节中仿真的一条股价路径绘于图形之上。

> S = c(100.0, 100.69, 99.74, 109.07, 103.68, 92.98, 101.40,          101.96, 114.84, 108.64, 108.63, 115.29, 110.33) > windows(width = 4.5, height = 3.3, pointsize = 8) > old.par <- par(mex = 0.8, mar = c(5,5,4,2) + 0.1) > plot(S, type="l")

 

我们先设定好股价向量S，一个有13个元素的向量，代表仿真一年走12步的路径。期初股价为100.0，期末价格为110.33。

下一行指令windows()表示设定绘图窗口宽为4.5英吋，高为3.3英吋，文字与符号字体大小为8。第三行指令是将旧的边界相关参数储存起来，并使用新的设定。mex=0.8是代表边界文字的缩放比，mar=c(5,5,4,2)+0.1表示下、左、上、右四个边界各留下5.1、5.1、4.1、2.1条线宽的距离。第四条指令则是将S的数据绘出。

 

一、常用参数

在plot指令之中，可以加入一些参数，改变绘图的输出效果。

xlim:x轴的范围，两元素的向量。

ylim:y轴的范围，两元素的向量。

xlab:x轴的名称。

ylab:y轴的名称。

main:图形的主标题。

sub:图形的副标题。

pch:点的型式，0-18为一般常用型式。

col:绘图的颜色。

 

Plot(S, xlim=c(0, 12), y=c(90, 120),xlab=”Step”, ylab=”Stock Price”,

main=”Monte-Carlo Simulation”, sub=”One Path Stock Price”, cex=0.8,

pch=16, col=”red”)

 

二、图形样式

在plot指令之中，type参数可以选择图形的样式。

type= “p” 点状图形

type= “l” 线状图形

type= “b” 点、线状不重迭图形

type= “o” 点、线状重迭图形

type= “h” 直线状分配图形

type= “n” 没有绘之空图，未来再增加

 

三、数学函数图

针对数学函数，可以使用curve()指令，传入函数名称、起始值、终止值，便可绘出图形。

> windows(width=4.5, height=3.3, pointsize=8)

> old.par <- par(mex=0.8, mar=c(5,5,4,2)+0.1)

> f <- function(x) x^2 - 3*x + 4 * sin(x)

> curve(f, from = -5, to = 5)