本文是课程《数据科学与金融计算》第4-5章的学习笔记，主要介绍金融数据处理、收益率计算和R与C++调用，用于知识点总结和代码练习，Q&A为问题及解决方案。
往期回顾：
| 博文 | 内容 |
|---|---|
| 【R】【课程笔记】01 R软件基础知识 | 数据类型、数据结构、运算、绘图等 |
| 【R】【课程笔记】02+03 基于R软件的计算 | 聚类分析、因子分析、神经网络、支持向量机等 |
| 【R】【课程笔记】04+05 数据预处理+收益率计算 | 金融数据处理、收益率、R与C++等 |
| 【R】【课程笔记】06 金融波动模型 | GARCH、SV、高频波动模型等 |
| 【R】【课程笔记】07 分位数回归与VaR（ES）计算 | VaR、ES、极值模型等 |
| 【R】【课程笔记】08 金融投资组合决策分析 | 均值-方差模型、均值-VaR模型、均值-CVaR模型等 |
目录
一、框架
二、代码
第四章 金融数据整理与预处理
4.2 金融数据格式
1、网页格式：xml、html，用“xml包”
2、从其他统计软件导入：命令“read.spss()”读取“.sav”格式SPSS文件，用“foreign包”
3、数据库格式：dbf，用“RODBC包”
4.3 金融数据导入
1、先进行变量设计和定义。
2、向量输入、数组输入、数据框输入（属性可以不一样，长度一样）、列表输入（每一个单独列出，$）
# 向量输入
(x<-c(11,12,13,14,15))
x<-scan()
# 数组输入
x<-1:12 #先给元素，再定义
a<-array(data=x,dim=c(3,4))
# 数据框输入
time<-c('2014-2-24','2014-2-25','2014-2-26','2014-2-27','2014-2-28')
AAPL.Volume<-c(10318200,8284000,9864900,10781500,13284600)
AAPL.Adjusted<-c(527.55,522.06,517.35,527.67,526.24)
AAPL<-data.frame(time,AAPL.Volume,AAPL.Adjusted)
# 列表输入
time<-c('2014-2-24','2014-2-25','2014-2-26','2014-2-27','2014-2-28')
AAPL.Volume<-c(10318200,8284000,9864900,10781500,13284600)
AAPL.Adjusted<-c(527.55,522.06,517.35,527.67,526.24)
mylist<-list(time=time,volume=AAPL.Volume,adjusted=AAPL.Adjusted)
print(mylist)
3、读取文件
# 加路径读取文件，不易出错
data <- read.csv('C:\\documents\\data.csv')
# 读取剪贴板
data <- read.delim('clipboard')
4.4 金融数据的预处理
1、时间序列的预处理
包括：数据合并（merge函数）、子集选择（subset(数据集，条件)）、随机抽样、数据补全（align函数）、频率转换、滚动窗操作（EWMA：指数滑动平均）。
library(timeSeries) # 时间序列
library(tseries) # 时间序列
library(RODBC) # 数据连接
library(quantmod) # 金融数量包
# (1) 数据合并
stock.names <- c('AAPL','GOOG')
getSymbols(stock.names, from='2014-03-01',to='2014-05-31') # 可能某些交易日没数据
start(AAPL) # 从哪开始
end(AAPL)
class(AAPL)
dim(AAPL) # 查看数据形状
dim(GOOG)
# 合并时间序列的函数，all=FALSE表示把两支股票都有的数据保留，TRUE表示所有数据都放进去（用NA填充）
assets <- merge.xts(AAPL, GOOG, all=FALSE)
dim(assets)
#保存R格式数据格式，方便下一次数据处理
save(assets, file='assets.RData')
load('assets.RData')
# (2) 子集选择
assets[1,]
assets <- assets[,c('AAPL.Close', 'GOOG.Close')] # 选择列
assets['2014-05'] # 选择行
assets[start(assets),] # 第一个记录
assets[end(assets),] # 最后一个记录
subset(assets, (assets[,1]>80) & (assets[,2]<550)) # 选择条件子集
# (3) 随机抽样
assets.ts <- as.timeSeries(assets) # 改变对象的类别，时间序列
class(assets.ts) # 查看类型
assets.samp <- sample(assets.ts, size=40) # 随机选样本，不放回地取
dim(assets.samp)
print(assets.samp)
sort(assets.samp) # 按时间排序
（1）数据补全：align函数
函数形式：align(x, by = “1d”, offset = “0s”,method = c(“before”, “after”, “interp”, “fillNA”,“fmm”, “periodic”, “natural”, “monoH.FC”),include.weekends = FALSE, …)
解释：interp：均值插值法
periodic：周期性，比如上个月
moniH.FC：单调法
# (4) 数据补全（随机抽取可能某些没有数据）
assets.ali <- align(assets.samp, by='1d', method='before', include.weekends=FALSE) #用前一天数据去补
dim(assets.ali)
print(assets.ali)
# (5) 频率转换
assets.m <- to.monthly(assets) #天数据转换为月数据
print(assets.m)
（2）滑动平均
①ans：对时间序列（applySeries），只需关注x、from、to三个属性。
②指数移动平均（EMA或EWMA）：指以指数式递减加权的移动平均。
# (6) 滑动窗
#applySeries对时间序列重复去做,按月求和
rollapply <- function(x, by, FUN, ...){
ans <- applySeries(x, from=by$from, to=by$to, by=NULL, FUN=FUN, format=x@format,
zone=finCenter(x), FinCenter=finCenter(x), title=x@title, documentation=x@documentation, ...)
attr(ans, 'by') <- data.frame(from=format(by$from), to=format(by$to))
ans
}
rts <- diff(log(assets)) # 直接对数差分收益率
rts <- returns(assets, method='continuous', percentage=TRUE) # 先填补再差分计算收益率
rts <- as.timeSeries(na.omit(rts)) # 删除MA，再变成时间序列
by <- periods(time(rts), period='1m', by='1d') # period:滑动的周期月（窗宽是1个月），by:按天记录的
rts.roll <- rollapply(rts, by=by, FUN='colSums') # colSums：按列和，收益率求和
print(rts.roll)
2、⭐⭐截面数据的预处理
包括：命名、标记数据、剔除数据、变量计算、数据合并和剔除NA值。
案例：要求剔除行业类型为金融保险行业的上市企业，剔除了样本区间所有性质在国有企业恶化其他类型企业之间
有变化的，剔除缺失数据的上市企业，剔除存在极端数据的样本观测。
Z <- odbcConnectExcel2007('FirmData.xlsx') # 数据比较大，先链接进来，运行更快
Data_Qvalue <- sqlFetch(Z, 'Qvalue', max=65536) # 直接读表格会比较慢
Data_Industry <- sqlFetch(Z, 'Industry', max=2717)
Data_Concent <- sqlFetch(Z, 'Concent', max=43227)
Data_Ownership <- sqlFetch(Z, 'Ownership', max=40899)
Data_Finance <- sqlFetch(Z, 'Finance', max=65532)
Data_Banlance <- sqlFetch(Z, 'Banlance', max=115619)
close(Z)
detach(package:RODBC)
tail(Data_Industry) # 检查一下
tail(Data_Concent)
tail(Data_Ownership)
tail(Data_Finance)
tail(Data_Banlance)
（1）股票命名方法：paste、substr
substr：取子串，从第“nchar(name.rev)-5”位取到第“nchar(name.rev)”位，最长6位。
# 命名唯一，用股票代码
firmName <- function(name){
name.rev <- paste('000000', as.character(name), sep='')
name.rev <- substr(name.rev, start=nchar(name.rev)-5, stop=nchar(name.rev))
name.c <- paste('c', name.rev, sep='')
name.c
}
Data_Qvalue$Comcd <- firmName(Data_Qvalue$Stkcd)
Data_Industry$Comcd <- firmName(Data_Industry$Stkcd)
Data_Concent$Comcd <- firmName(Data_Concent$Stkcd)
Data_Ownership$Comcd <- firmName(Data_Ownership$Stkcd)
Data_Finance$Comcd <- firmName(Data_Finance$Stkcd)
Data_Banlance$Comcd <- firmName(Data_Banlance$Stkcd)
（2）删减数据：判断
A%in%B：表示如果A里的元素在B里，在则返回TRUE
# 去掉金融行业
Data_Industry <- Data_Industry[!is.na(Data_Industry$Indcd),] #去掉NA
#不等于金融行业（1）的拿出来
NonfinComcd <- levels(factor(as.character(Data_Industry[Data_Industry$Indcd != 1, 'Comcd'])))
# %in%表示如果前面的在后面就取出来：即：把非金融的行业拿出来
Qvalue.Nonfin <- Data_Qvalue[(Data_Qvalue$Comcd %in% NonfinComcd), ]
Industry.Nonfin <- Data_Industry[(Data_Industry$Comcd %in% NonfinComcd), ]
Concent.Nonfin <- Data_Concent[(Data_Concent$Comcd %in% NonfinComcd), ]
Ownership.Nonfin <- Data_Ownership[(Data_Ownership$Comcd %in% NonfinComcd), ]
Finance.Nonfin <- Data_Finance[(Data_Finance$Comcd %in% NonfinComcd), ]
Banlance.Nonfin <- Data_Banlance[(Data_Banlance$Comcd %in% NonfinComcd), ]
# 重新命名
names.Qvalue <- c('Comcd', 'Accper', 'QVal')
names.Industry <- c('Comcd', 'Indcd')
names.Concent <- c('Comcd', 'Reptdt', 'OwnCon1', 'OwnCon5', 'OwnCon10')
names.Ownership <- c('Comcd', 'Reptdt', 'Indicator', 'State')
names.Finance <- c('Comcd', 'Accper', 'ROA', 'ROE')
names.Banlance <- c('Comcd', 'Accper', 'CurrentAsset', 'TotalAsset', 'CurrentLib', 'TotalLib')
#取出第6-7个字符（即月份），即取出月份是12的数据，拿出年终数据
Qvalue.tab <- Qvalue.Nonfin[substr(Qvalue.Nonfin$Accper, start=6, stop=7)=='12', names.Qvalue]
Industry.tab <- Industry.Nonfin[,names.Industry]
Concent.tab <- Concent.Nonfin[substr(Concent.Nonfin$Reptdt, start=6, stop=7)=='12', names.Concent]
Finance.tab <- Finance.Nonfin[substr(Finance.Nonfin$Accper, start=6, stop=7)=='12', names.Finance]
Banlance.tab <- Banlance.Nonfin[substr(Banlance.Nonfin$Accper, start=6, stop=7)=='12', names.Banlance]
#所有制类型
Ownership.tab <- Ownership.Nonfin[Ownership.Nonfin$Indicator==1, names.Ownership]
Ownership.tab$state <- 0
Ownership.tab$state[(Ownership.tab$State == 1100) | (Ownership.tab$State == 2100)] <- 1
（3）变量计算
#用Accper的第1-4字符赋值给Enddt
Qvalue.tab$Enddt <- as.numeric(substr(Qvalue.tab$Accper, 1, 4))
Concent.tab$Enddt <- as.numeric(substr(Concent.tab$Reptdt, 1, 4))
Ownership.tab$Enddt <- as.numeric(substr(Ownership.tab$Reptdt, 1, 4))
Finance.tab$Enddt <- as.numeric(substr(Finance.tab$Accper, 1, 4))
Banlance.tab$Enddt <- as.numeric(substr(Banlance.tab$Accper, 1, 4))
#股权制衡度（EBD）、企业规模、流动比例
Concent.tab$EBD <- Concent.tab$OwnCon5/Concent.tab$OwnCon1
Banlance.tab$Size <- log(Banlance.tab$TotalAsset)
Banlance.tab$Currt <- Banlance.tab$CurrentAsset/Banlance.tab$CurrentLib
Banlance.tab$AssLibRatio <- Banlance.tab$TotalLib/Banlance.tab$TotalAsset
（4）合并
KeyField1 <- 'Comcd' # 按照公司名称
KeyField2 <- 'Enddt' # 按照年份
#两两去做合并
Qvalue.Finance <- merge(Qvalue.tab, Finance.tab, by=c(KeyField1, KeyField2))
Qvalue.Finance.Concent <- merge(Qvalue.Finance, Concent.tab, by=c(KeyField1, KeyField2))
Qvalue.Finance.Concent.Ownership <- merge(Qvalue.Finance.Concent, Ownership.tab, by=c(KeyField1, KeyField2))
Qvalue.Finance.Concent.Ownership.Banlance <- merge(Qvalue.Finance.Concent.Ownership, Banlance.tab, by=c(KeyField1, KeyField2))
Qvalue.Finance.Concent.Ownership.Banlance.Industry <- merge(Qvalue.Finance.Concent.Ownership.Banlance, Industry.tab,
by=c(KeyField1))
# 移除NA值
var.names <- c("Comcd", "Enddt", "QVal", "ROA", "ROE", "OwnCon1", "OwnCon10", "EBD", "state","Size",
"CurrentAsset", "TotalAsset", "CurrentLib", "TotalLib", "Currt", "AssLibRatio", "Indcd")
Data <- Qvalue.Finance.Concent.Ownership.Banlance.Industry[ ,var.names]
Data <- na.omit(Data)
# 设置年份的虚拟变量
Data$Year02 <- 0
Data$Year03 <- 0
Data$Year04 <- 0
Data$Year05 <- 0
Data$Year06 <- 0
Data$Year07 <- 0
Data$Year08 <- 0
Data$Year09 <- 0
Data$Year10 <- 0
Data$Year11 <- 0
Data$Year12 <- 0
Data$Year13 <- 0
Data[Data$Enddt == 2004, 'Year04'] <- 1
Data[Data$Enddt == 2005, 'Year05'] <- 1
Data[Data$Enddt == 2006, 'Year06'] <- 1
Data[Data$Enddt == 2007, 'Year07'] <- 1
Data[Data$Enddt == 2008, 'Year08'] <- 1
Data[Data$Enddt == 2009, 'Year09'] <- 1
Data[Data$Enddt == 2010, 'Year10'] <- 1
Data[Data$Enddt == 2011, 'Year11'] <- 1
Data[Data$Enddt == 2012, 'Year12'] <- 1
Data[Data$Enddt == 2013, 'Year13'] <- 1
dim(Data)
print(Data)
第五章 金融资产收益计算
多期简单收益率（连乘）：
多期连续复利收益率（求和）：
红利收益率、超额收益率
5.2 股票类资产收益率
单个股票收益率：算术平均值、几何平均值
（1）线性关系、乘积关系
（2）连续复利、简单收益率
prod：连乘
library(zoo) # 时间序列
library(timeSeries)
# 1. 计算收益
# (1) 读数据
da <- read.table("5-1.txt",head=T,colClasses="character")
clsprc <- as.numeric(da[,3])
dates <- as.Date(da[,2])
cls <- zoo(clsprc, dates) # 改变对象类型
lclsprc <- lag(cls,k=-1) # 延迟一期，-1删除最大的那个
clsprc <- as.numeric(da[2:237,3])
lclsprc <- as.numeric(lclsprc)
# (2) 计算日收益
re <- 100 * (log(clsprc)-log(lclsprc))
re.1 <- 100 * diff(log(cls))
re.2 <- returns(cls, method='continuous', percentage=TRUE)
re[1:10]
re.1[1:10]
re.2[1:11,]
# (3) 算术平均
(am <- mean(re))
(am.1 <- mean(re.1))
(am.2 <- mean(re.2, na.rm=TRUE))
# (4) 几何平均
#(gm <- 100*(cumprod(re/100+1)[length(cumprod(re/100+1))]^(1/length(re))-1))
100*((prod(re/100+1))^(1/length(re))-1)
多个股票收益率计算
da1 <- read.table("CAQC.txt",head=T,colClasses="character")
da2 <- read.table("ZXTX.txt",head=T,colClasses="character")
wcls1 <- as.numeric(da1[,3])
wcls2 <- as.numeric(da2[,3])
dates1 <- as.yearmon(da1[,2])
dates2 <- as.yearmon(da2[,2])
cls1 <- zoo(wcls1,dates1)
cls2 <- zoo(wcls2,dates2)
dat.merge <- merge(cls1, cls2)
# 计算收益率
r.merge <- returns(dat.merge, method='continuous', percentage=TRUE)
print(r.merge)
资产组合收益率计算：关键在比例
·
# 3. calculate portfolio returns
# (1) read data
da <- read.table("5-3.txt",head=T)
sum(is.na(da))
p.app <- na.approx(da)
r.stocks <- diff(log(p.app))*100
# (2) 模拟组合权重
set.seed(12345)
w <- runif(ncol(da), 0, 1)
(w <- w/sum(w))
# (3) 投资组合收益
rp <- r.stocks %*% w
5.3 债券类资产收益率
三种收益计算：
1、到期收益率：uniroot(函数，解的范围c(0,1)，传的参数p，Cs)
# (1) 定义函数
f <- function(r, p, Cs){
n <- length(Cs)
tt <- 1:n
loss <- p - sum(Cs/((1+r)^tt))
loss
}
# (2) 求解
Cs <- c(第一期现金流, 第二期, 第三期, 第四期)
p <- 初始价格
uniroot(f, c(0,1), p=p, Cs=Cs)
含本金：修改函数
f_exp <- function(r, p, Cs, Cp){
n <- length(Cs)
tt <- 1:n
loss <- p - sum(Cs/((1+r)^tt)) - Cp/(1+r)^n
loss
}
Cs <- rep(200, times=40)
Cp <- 2500
p <- 2364
uniroot(f_exp, c(0,1), p=p, Cs=Cs, Cp=Cp)
比较期望收益率与内生收益率
# 函数同上
Cs <- rep(950000, times=15)
Cp <- 15000000
p <- 22710000
r.expected <- 0.067
r.half <- uniroot(f_exp, c(0,1), p=p, Cs=Cs, Cp=Cp)$root #取列表
r.annulized <- r.half*2 #年化
if (r.annulized>=r.expected){
cat('Do it, and the endougeneous return is', r.annulized, '\n')
}else{
cat('Not do it, and the endougeneous return is', r.annulized, '\n')
}
2、有效年利率
r/m：周期性利率
# (1) 定义函数
eff_rts <- function(pr,m){
er<-(1+pr)^m-1
}
# (2) 计算
(da1 <- eff_rts(0.05,2))
(da2 <- eff_rts(0.025,4))
3、债券久期与凸度
久期度量的是投资回收的平均时间，而修正久期实质上度量的是债券价格对收益率的一阶导数[1]。
（1）修正久期：
修正久期大的债券，利率上升所引起其价格下降幅度就越大，而利率下降所引起其价格上升幅度也越大。因此，修正久期大的债券比修正久期小的债券抗利率下降风险能力强，但抗利率上升风险能力较弱。
# （1）定义函数
jq1 <- function(y,coupon,period,p0){
c2<-0
tc2<-0
for(n in 1:period){
t=n
if(n<period) c<-coupon else if(n==period) c<-coupon+p0
c1=c/(1+y)^n
tc1<-t*c1
c2<-c2+c1
tc2<-tc2+tc1
}
d<-tc2/(c2*2) # duration
md<-d/(1+y) # adjusted duration
list(d=d,md=md)
}
# （1）修改后的函数
jq2 <- function(y,coupon,period,p0){
c2<-0
tc2<-0
coupon = rep(coupon, period)
coupon[period] = coupon[period]+p0
for(n in 1:period){
t1=n
c1=coupon[n]/(1+y)^n
tc1<-t1*c1
c2<-c2+c1
tc2<-tc2+tc1
}
d<-tc2/(c2*2) # duration
md<-d/(1+y) # adjusted duration
list(d=d,md=md)
}
# (2) 计算
(re1<-jq1(0.05,5,10,2000))
(re2<-jq2(0.05,5,10,2000))
（2）凸度：债券价格对收益率的二阶导（变化速度）
# (1) 定义函数
td <- function(y, cupon, period, p0){
c2<-0
tc2<-0
for(n in 1:period){
t=n
if(n<period) c<-cupon else if(n==period) c<-cupon+p0
c1=c/(1+y)^n
tc1<-t*(t+1)*c1
c2<-c2+c1
tc2<-tc2+tc1
}
covex<-tc2/(c2*((1+y)^2))
yearlycovex<-covex/4
list(covex=covex,yearlycovex=yearlycovex)
}
# (2) do calculation
(re <- td(0.05,5,10,1500))
# (3) calculate the impact of convexity
delta.r <- 0.04
(effects <- 0.5*re$yearlycovex*delta.r^2)
5.4 收益率的分布及其特征
1、统计分布：联合分布、边缘分布、条件分布
2、数字特征（函数：basicStats）：样本均值、样本方差、样本偏度、样本峰度、分位数
3、常用分布函数：正态分布、对数正态分布、混合正态分布、稳态分布//
重尾分布[2]：广义极值分布
4、多元收益率统计
似然函数、多元正态分布、条件密度函数
library(fBasics) # for basic statistics
library(mnormt) # for multivariate normal distribution
library(ks) # for kernel smooth
da <- read.table('5-21.txt',head=T)
dim(da)
rt <- da[,2:3]
basicStats(da[,2]) # 描述性统计量
# 模拟多元正态分布
(m1 <- apply(rt,2,mean))
(v1 <- cov(rt))
set.seed(12345)
x <- rmnorm(1000, mean=m1,varcov=v1)
dim(x)
# 比较
par(mfrow=c(1,2))
plot(rt[,2], rt[,1], xlab='上证', ylab='深证', cex=0.8, main='真实数据')
abline(a=0, b=1, lty=2)
plot(x[,2], x[,1], xlab='上证', ylab='深证', cex=0.9, main='模拟数据')
abline(a=0, b=1, lty=2)
par(mfrow=c(1,1))
# (4) estimate kernel density of real data
dens.hat <- kde(rt)
plot(dens.hat, display='filled.contour2', cont=seq(10,90,by=10), xlim=c(-0.06, 0.06), ylim=c(-0.05, 0.06))
plot(dens.hat, display='persp',zlab='密度函数', theta=30, phi=60,thin=10, border=1, col="white")
补充：R与C++结合
1、Gibbs抽样
（1）C++代码：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix gibbs_cpp(int N, int thin) {
NumericMatrix mat(N, 2);
double x = 0, y = 0;
for(int i = 0; i < N; i++) {
for(int j = 0; j < thin; j++) {
x = rgamma(1, 3, 1 / (y * y + 4))[0];
y = rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))[0];
}
mat(i, 0) = x;
mat(i, 1) = y;
}
return(mat);
}
（2）R代码：
gibbs_r <- function(N, thin) {
mat <- matrix(nrow = N, ncol = 2)
x <- 0
y <- 1
for (i in 1:N) {
for (j in 1:thin) {
y <- rexp(1)
x <- rnorm(1, mean(x), y)
}
mat[i, ] <- c(x, y)
}
mat
}
2、循环：f(n)=f(n-1)+f(n-2)
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
int fib_cpp_1(int n)
{
if(n==1||n==2) return 1;
return fib_cpp_1(n-1)+fib_cpp_1(n-2);
}
// [[Rcpp::export]]
double sumC(NumericVector x) {
int n = x.size();
double total = 0;
for(int i = 0; i < n; ++i) {
total += x[i];
}
return total;
}
// [[Rcpp::export]]
NumericVector rowSumsC(NumericMatrix x) {
int nrow = x.nrow(), ncol = x.ncol();
NumericVector out(nrow);
for (int i = 0; i < nrow; i++) {
double total = 0;
for (int j = 0; j < ncol; j++) {
total += x(i, j);
}
out[i] = total;
}
return out;
}
r代码
cppFunction(
'int fib_cpp_0(int n){
if(n==1||n==2) return 1;
return(fib_cpp_0(n-1)+fib_cpp_0(n-2));
}'
)
fib_r <- function(n){
if(n==1||n==2) return(1)
return(fib_r(n-1)+fib_r(n-2))
}
sourceCpp("fib.cpp")
⭐⭐例题：随机抽样
（1）C++代码：
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix random_cpp(int N, int M) {
NumericMatrix mat(N, 2);
mat(0,0)=0;
mat(0,1)=1;
double x=0, y=0;
for(int i=1; i<N; i++){
for(int j=0; j<M; j++){
y=rexp(1, 1)[0];
double x_sum=0;
for(int t=0; t<i; t++){
x_sum+=mat(t,0);
}
x=rnorm(1, x_sum/i, sqrt(y))[0];
}
mat(i,0)=x;
mat(i,1)=y;
}
return(mat);
}
（2）R代码 & 调用C++
library(Rcpp)
# 用R程序编写生成随机数
random_r<-function(N,M){
mat<-matrix(nrow=N,ncol=2)
mat[1,1]<-0
mat[1,2]<-1
x<-y<-0
for(i in 2:N){
for(j in 1:M){
y<-rexp(1, rate=1)
x<-rnorm(1, mean(mat[1:i-1,1]), sqrt(y))
}
mat[i,]<-c(x,y)
}
mat
}
#以N=100为例，查看R程序编写生成随机数
set.seed(123)
samp1=random_r(100,100)
dim(samp1)
samp1[1:20,]
#用R软件调用编写C++程序生成随机数
sourceCpp("rand.cpp") #调用C++程序"rand.cpp"计算随机数
set.seed(123) #为了比较R与C++生成的随机数是否相同，设置相同的种子数
samp2=random_cpp(100,100)
samp2[1:20,]
dim(samp2)
#比较两种方法的运行时间和结果#
set.seed(123)
system.time(random_r(1000,1000))
system.time(random_cpp(1000,1000))
2957
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
