朱世武《基于SAS系统的金融计算》书中功能用R的实现

abel

国内金融数据处理方面用SAS的恐怕朱世武教授的两本书是听过的，初次见到他的书的时候，就动过念头——将其中的代码换成R版本的（金融计算用C++、Matlab、Excel VBA等的很多，前者是“高级”的，后者是“大众化”的选择）。前年的时候，还计划把这想法结合教学实践来展开下。不过，自己太懒了，落下了。现在开个帖子，争取每周都小更新下，大家督促，并且修改其中的错误，进行优化。

这两书名分别为：《SAS编程技术与金融数据处理》（2003）《基于SAS系统的金融计算》（2004）。前者简要介绍了一些金融数据库的内容，主题部分是SAS编程的一些基本技术（包含了简单的MARCO和IML的介绍），该书最后两章涉及了金融计算的入门技术——资产收益率、固定收入证券的计算。后书的前面三章总体上和前书一样，后面的章节逐次介绍了指数编制、波动率计算等金融计算的各个方面。

就我个人的角度来说，我不太习惯SAS做金融计算，感觉不够灵活，除非针对整理好了的做固定的模型，否则我不会选择SAS作为干活的工具。如果看过朱教授的这两本书以后，估计你也会是这种感觉——甚至有人会愿意去学C++去了。相对来说，R是一个非常不错的选择（这里并不涉及R在CRAN上众多的pkgs，而只适用其二三十个左右推荐包的范围；因为R中已经有若干来做金融计算和建模的pkgs）：容易上手、语法简便、原型开发以后可以迅速的转向C++，etc。（我是R的粉丝，遇到有人计算比较犯难的时候，就推荐其尝试下R）。

abel

NOTE：2003这本书从第十九章开始（该书前面介绍了金融数据库和纯SAS技术，所以暂时就不给出了）

2004这本书的第二章和这章内容一致。

首先，提供Example19.2的另外一种实现方法（在SAS中对不同record的数据进行计算很麻烦，比较好的选择是使用retain，下面就是这个思路）

<br />
data aveRtn(drop=prod sum index);<br />
    set a1;<br />
    retain prod 1 sum 0 index 0;<br />
    index = index + 1;<br />
    sum=sum + r_log;<br />
    prod = prod * (1 + r_pct);<br />
    result_pct=prod**(1/index) - 1;<br />
    result_log= sum/index;<br />
run;<br />
<br />
<br />
# Example 19.1<br />
dat <- read.csv("a1a0001.csv", stringsAsFactors=FALSE)[, c("date", "Clpr")]; head(dat)<br />
dat$date <- as.Date(dat$date);<br />
dat <- dat[order(dat$date),]<br />
<br />
# 以1990年为基年计算月份和季度，如计算week的方式，均可转化为julian历法后计算<br />
year <- as.integer(substr(dat$date, 1, 4))<br />
qtr <- (year-1990)*4 + 1 + floor((as.integer(substr(dat$date, 6, 7))-1)/3)<br />
month <- (year-1990)*12 + as.integer(substr(dat$date, 6, 7))<br />
week <- as.integer(3 +julian(dat$date) - julian(dat$date[1]))%/%7    # 其实此法最简便！<br />
<br />
dat.year <-dat[cumsum(tapply(year, year, length)), ]<br />
dat.qtr <- dat[cumsum(tapply(qtr, qtr, length)), ]<br />
dat.month <- dat[cumsum(tapply(month, month, length)), ]<br />
dat.week <- dat[cumsum(tapply(week, week, length)), ]<br />
<br />
rtn <- function(price){<br />
    x <- price<br />
    r_pct <- c(NA, diff(x)/x[-length(x)])<br />
    r_log <- c(NA, diff(log(x)))<br />
    ans <- data.frame(r_pct, r_log)<br />
    ans<br />
}<br />
ret.year<- cbind(dat.year, rtn(dat.year$Clpr))<br />
ret.qtr<- cbind(dat.qtr, rtn(dat.qtr$Clpr))<br />
ret.month<- cbind(dat.month, rtn(dat.month$Clpr))<br />
ret.week <- cbind(dat.week, rtn(dat.week$Clpr))<br />
<br />
# Example 19.2<br />
ret.mean <- function(dat){<br />
    dat <- dat<br />
    year <- as.integer(substr(dat$date, 1, 4))<br />
    ret.pct <- dat[year>= 1995 & year <=2000, "r_pct"]<br />
    gm <- (prod(1+ret.pct))^(1/6) - 1<br />
    am <- mean(ret.pct)<br />
    c(gm, am)    <br />
}<br />
ret.mean(ret.year)<br />
ret.mean(ret.qtr)<br />
ret.mean(ret.month)<br />
ret.mean(ret.week)<br />
<br />
<br />
# Example 19.3 <br />
#    Clpr        收盘价|Close Price<br />
#    Rigrat    配股比例|right shares ratio<br />
#    date        日期|Date<br />
#    divamt    每股派息额|Dividend Cash Amount<br />
#    divrat    每股送股比例| dividends ratio<br />
#    reispr    增发新股股价|reissue price<br />
#    reisvol    增发新股数量|Reissue volume<br />
#    rigobv    实际配股数|right shares observed<br />
#    rigpr        配股价|Right Shares Price<br />
#    shrout    观测到的总股本|Shares Outstanding<br />
dat <- read.csv("a600601.csv", stringsAsFactors=FALSE)[, c("date", "Clpr", "clpr_r", <br />
    "divrat", "Rigrat", "reisvol", "shrout", "rigpr", "reispr", "divamt")]; head(dat)<br />
index <- read.csv("a1a0001.csv", stringsAsFactors=FALSE)[, c("date", "sums")]; <br />
dat <- merge(dat, index, all.y=TRUE)    # 个股停牌处理<br />
year <- as.integer(substr(dat$date, 1, 4))<br />
dat <- dat[year>=1997 & year<=2000, ]<br />
<br />
attach(dat)<br />
rtn.daily <- (rtn(Clpr))[, 2]<br />
inds <- apply(dat[, c(4:6, 8:10)], 1, function(x)!any(is.na(x)))<br />
inds <- (1:nrow(dat))[inds]<br />
for(ind in inds){<br />
    clpr_adj <- Clpr[ind]*(1 + divrat[ind] + Rigrat[ind] + reisvol[ind]/shrout[ind-1]) <br />
            - rigpr[ind]*Rigrat[ind] - reispr[ind]*reisvol[ind]/shrout[ind-1] + divamt[ind]<br />
    rtn.daily[ind] <- log(clpr_adj) - log(Clpr[ind-1])<br />
}<br />
dat$r.daily <- rtn.daily<br />
detach(dat); remove(shrout)<br />
dat <- dat[, c("date", "Clpr", "clpr_r", "r.daily")]<br />
<br />
# Example 19.4 ARIMA(0, 1, 0)<br />
set.seed(123456)<br />
eps <- c(0, rnorm(1050, 0, 1))<br />
p.t <- 8.6 + cumsum(eps)<br />
p.t <- p.t[51:1050]<br />
ts.plot(p.t)<br />
<br />
# Example 19.5 ARIMA(0, 1, 0)<br />
set.seed(123456)<br />
eps <- c(0, rnorm(1050, 0, 1))<br />
r.t <- -0.02 + 0.5*cumsum(eps)<br />
r.t <- r.t[51:1050]<br />
ts.plot(r.t)<br />
<br />
# Example 19.6 <br />
set.seed(123456)<br />
mu <- 0<br />
eps <- c(0, rnorm(1050, 0, c(rep(1, 110), rep(2, 200))))<br />
r.t <- -0.02 + mu + 0.5*cumsum(eps)<br />
r.t <- r.t[11:310]<br />
ts.plot(r.t)<br />
<br />
# Example 19.7 ARIMA(0, 1, 1)<br />
set.seed(123456)<br />
eps <- rnorm(1050)<br />
eps <- c(NA, eps[2:length(eps)] - 0.8 * eps[-length(eps)])<br />
r.t <- -0.01 + cumsum(eps[-1])<br />
r.t <- r.t[51:1050]<br />
ts.plot(r.t)<br />
# 其实都可以看成是ARIMA模型，直接用arima.sim函数就能搞定，如此编写是为了和书中保持相同思路<br />
<br />
# Example 19.8 GARCH(1, 1)<br />
t <- 1:1052<br />
e <- c(0, rnorm(1050))<br />
eps <- h<- r.t <- rep(0, 1051)<br />
for(i in 2:1052){<br />
    h[i+1] <- 0.08 + 0.3 * (e[i])^2 + 0.3 * h[i]<br />
    eps[i+1] <- sqrt(h[i+1]) * e[i+1]<br />
    r.t[i+1] <- eps[i+1] + 0.5 * r.t[i] - 0.03 * r.t[i-1]<br />
}<br />
r.t <- r.t[51:1050]<br />
ts.plot(r.t)<br />
<br />
# The END<br />

未完待续 ……

abel

2003书中第二十章和2004书中第三章内容一致。

本章包含固定收入证券计算的一些内容，有固定收入证券的收益率、债券定价、债券的必要收益率等内容，此外还包含久期、凸度计算，也有二叉树定价模型。

在计算收益率的时候，核心的问题是要处理好一个高次方程解的问题，就是一个现金流量在该收益率的情况下现值为0了，这个高次方程可以从数学上证明通常都可以通过迭代的方法来求解。

在朱教授的书中，使用了所谓的”试错法“，实际上就是”人肉“逼近法。其实在解决该类问题已经有很多好用的自动迭代的方法没有必要”人肉“。比如二分法、线性插值法、Newton-Raphson法等，下面代码给出了牛顿法迭代求解（该方法性质那是相当的好），以及使用R中已有的polyroot函数——请注意，在持续期数太多的时候，解这个高次方程无疑是很恐怖的事情！不相信的同志可以试一试是的一个期数超过1000的方程玩玩，保管印象会深刻些。

<br />
# Newton-Raphson method, used in Example 1 to 5<br />
irr.NR <- function(cf, init=1, itor=20){<br />
    if(cf[1]>0) cf[1] <- -cf[1]    # Cash out at first<br />
    x <- rep(init, itor)<br />
    k <- length(cf) - 1<br />
    for(i in 1:itor){    # Newton-Raphson method<br />
        x[i+1] <- x[i] - sum(cf * (x[i])^(0:k))/(sum(cf[-1] * (1:k) * (x[i])^(0:(k-1))))<br />
    }<br />
    return(1/x[itor] - 1)<br />
}<br />
<br />
# polyroot method, used in Example 1 to 5<br />
irr.poly <- function(cf){<br />
    if(cf[1]>0) cf[1] <- -cf[1]    # Cash out in the point!<br />
    roots <- polyroot(cf)<br />
    roots <- Re(roots[abs(Im(roots))<1e-7])<br />
    roots <- 1/roots - 1<br />
    return(sort(roots[roots>=0 & roots<=1]))<br />
}<br />
<br />
# Example 20.1<br />
cf <- c(7704, 2000, 2000, 2500, 4000)<br />
irr.NR(cf)<br />
irr.poly(cf)<br />
<br />
# Example 20.2<br />
cf <- c(1243.82, rep(50, 9), 1050)<br />
irr.NR(cf)<br />
irr.poly(cf)<br />
<br />
# Example 20.2<br />
cf <- c(944.14, rep(50, 9), 1050)<br />
irr.NR(cf)<br />
irr.poly(cf)<br />
<br />
# Example 20.4, Note: when the time last too long, the polyroot consumes so much time!<br />
cf <- c(43449, rep(349.6, 360))<br />
irr.NR(cf)<br />
irr.poly(cf)<br />
<br />
# Example 20.5<br />
cf <- c(19696024, rep(1000000, 29), 21000000)<br />
irr.NR(cf)<br />
irr.poly(cf)<br />

紧接着我们就要解决真实利率（有效利率，effective rate）和名义利率(nominal rate)之间相互转换的关系，这里顺带也给出了利息强度（连续复利）的转换。同时提醒一下同志们，原书中计算债券的那个利率的描述其实是有点小问题的——当然这里不纠结利息理论的那些东西了，因为从数学上来说是忒简明的东西。

<br />
# effective and nominal rates transformation methods<br />
rate.eff <- function (nom, p)    # non 2 eff<br />
{<br />
    if (p != Inf) {    # Force of interest<br />
        return((1 + nom/p)^p - 1)<br />
    }<br />
    else {<br />
        return(exp(nom) - 1)<br />
    }<br />
}<br />
rate.nom <- function (eff, p)        # eff 2 non<br />
{<br />
    if (p != Inf) {    # Force of interest<br />
        return(p*((1+eff)^(1/p)-1))<br />
    }<br />
    else {<br />
        return(log(eff + 1))<br />
    }<br />
}<br />
<br />
# Example 20.6<br />
rate.eff(0.04*2, 2)<br />
rate.eff(0.02*4, 4)<br />

其实编写下面这个债券收益率函数的时候，我纠结了好久，一般习惯吧债券搞成半年给一次票息，于是名称上就有点小混乱了。不过把债券理解为固定（一般为等时间间隔上）的收入流就万事大吉了。（以后很好奇为啥通常都是半年给，而不是一年或一个月给，至少年利率和月利率是比较容易理解的嘛）

<br />
bond.irr <- function(issue=1000, coupon=0.03, maturity=60, price=1000){<br />
    cf <- c(-price, rep(coupon*issue, maturity-1), (1+coupon)*issue)<br />
    return(2*irr.NR(cf))    # return(rate.eff(2*irr.NR(cf), 2)) -> effective rates<br />
}<br />
# Example 20.7, bond's irr<br />
bond.irr(issue=1000, coupon=0.03, maturity=36, price=700.89)<br />

下个跟帖将讨论下非完整年的度量问题，到时候顺便痛骂下无良的金融资本家一下，未完待续……

paladin_aj

Excellent!!!

waiting for the following contents.....

emilio

哇，好强，支持支持

danielbaggio

强烈支持，我也是R的爱好者！

朱世武的书是挺好，但我老觉得SAS不是用来搞金融的。。。当然SAS我连皮毛都不懂，个人感觉。。。

另外，朱世武第二本书里面VaR的back testing，不知道R有没有能实现的包？

期待楼主不断更新，我也常来关注楼主更新，嘿嘿

abel

固定收益证券（fixed-income instrument），顾名思义是在特定的时间内取得固定的收益并预先知道取得收益的数量和时间的一种证券。简单

理解，就是一组确定的现金流（每次付账的时间和额度是确定的）。所以从数学的角度来看，固定收益证券的价格收益率等应该是很容易计算

的——就是一个数学方程能表达的模型。当计息时间和付息时间一致的时候（或者是整倍数），比较容易理解和计算。当计息时间和付息时间（

清算日处于两个到期日之间）不一致的时候，就出现非整数的年，比如0.2145年之类的。

值得注意的是，我们现在使用的Julian历法中每年和每月的天数是不相等的！如果没闰年、大小月、闰月，我们就可以直接某某天数除以365，

就可以得到一个小数。但是因为有上面提到的麻烦在，计算非整数年的时候，就有点小纠结了。总体上来说，有严格计算实际天数、一年记为

固定的360、365天，一个月记为30天等，加上月末如何计算，规则还是比较多，不过通常的方法有如下几种：(来自http://en.wikipedia.org/wiki/Day_count_convention)

1、30/360法则

    1.1、30/360 US

    1.2、30E/360

    1.3、30E/360 ISDA

    1.4、30E+/360

2、实际法则

    2.1、实际/实际 ICMA

    2.2、实际/实际 ISDA

    2.3、实际/365 Fixed

    2.4、实际/360

    2.5、实际/365L

    2.6、实际/实际 AFB

其中的ISDA ICMA AFB等都是一些组织，计算天数的方法就要依照其规定来。通常我们遇到的有三类：30/360、实际/360、实际/实际三种，其中前两种，相对来说付利息的人就有点小吃亏了，大家有兴趣的可以自己计算下。（以前有新闻报道说我们在银行取利息的时候，就采用对客户不利的那种来，貌似大家要看清楚银行的规定吧）

那么我们就转入下面一个实际的例子：

票息率10%的公司债，2003年3月1日到期。该债券的价格为118.788元，清算日为1997年7月17日。要计算到期收益率。

1、非整年数=清算日到下一个票息日的天数除以票息天数=44/180=0.244444 （30/360法则）

2、解方程组 -118.788 + 50*v^(0.244) + 50*v^(1.244)+...+50*v(10.244)+ 105*v^(11.244)=0就得到了折现因子v，收益率=1/v - 1

现金流为：

cf <- cf <- c(-118.788, rep(5, 11), 105)

现金流发生的时间为：

para <- c(0, w + 0:11)

依据朱书上得到的结果：3.63%，有

cf * (1/1.0363)^para

-2.560512                # 注意，这里距离等于0还差比较多！

修改后的irr.NR函数

<br />
irr.NR.factor <- function(cf, init=1, itor=10, w=1){<br />
    if(cf[1]>0) cf[1] <- -cf[1]# Cash out in the point!<br />
    x <- rep(init, itor)<br />
    para <- c(1-w, 1:(length(cf)-1))<br />
    for(i in 1:itor){<br />
        x[i+1] <- x[i] - sum(cf * (x[i])^para)/sum(para * cf * x[i]^(para-1))<br />
    }<br />
    return(1/x[itor] - 1)<br />
}<br />
irr.NR.factor(cf, w=44/180)<br />
0.03372264<br />
# 验算：<br />
sum(cf * (1/1.03372264)^para)<br />
2.524802e-06        # 误差就相对小不少了<br />

注意：债券收益率计算的时候，市场价格有的是净价格，有的是肮脏价格，采用两种不同价格计算的结果不一，这里不涉及金融理论，仅作通过计算示例。

abel

投资组合的内生收益率计算

1、确定投资组合中的现金流

2、建立方程式，求解

Example 3.9

三种债券，票息支付日相同。投资组合市值为57 259 000元。每只债券的现金流量为：

A <- c(rep(350000, 9), 10350000)

B <- c(rep(1050000, 13), 21050000)

C <- c(rep(900000, 6), 30900000)

len <- 14

port <- cbind(c(A, rep(0, len-length(A))), c(B, rep(0, len-length(B))), c(C, rep(0, len-length(C))))

colnames(port) <- LETTERS[1:3]

port.cf <- apply(port, 1, sum)

cbind(port, port.cf)

             A        B        C  port.cf

[1,]   350000  1050000   900000  2300000

[2,]   350000  1050000   900000  2300000

[3,]   350000  1050000   900000  2300000

[4,]   350000  1050000   900000  2300000

[5,]   350000  1050000   900000  2300000

[6,]   350000  1050000   900000  2300000

[7,]   350000  1050000 30900000 32300000

[8,]   350000  1050000        0  1400000

[9,]   350000  1050000        0  1400000

[10,] 10350000  1050000        0 11400000

[11,]        0  1050000        0  1050000

[12,]        0  1050000        0  1050000

[13,]        0  1050000        0  1050000

[14,]        0 21050000        0 21050000

(irr <- irr.NR((cf <- c(-57259000, port.cf))))

结果为：0.0468378

验算：

sum(cf * (1/(1+irr))^(0:14)) # 3.993046e-08

而原书中结果为：0.0476966

sum(cf * (1/(1+0.0476966))^(0:14)) # 3.993046e-08

验算结果为： -383137.3，误差比较大

最好我们可以得到半年计息的年名义收益率为：0.0468378*2= 0.0936756

小结：固定收益证券计算收益率的时候，确定每个时点的现金流量，时点间的时间差，遇到不为整年的情况就依照该证券的规定来处理，这样就可以得到一个方程，使用迭代法求解该房产近似值，然后验算。迭代法推荐使用牛顿法（只在两类情况下，该方法会歇菜，但是上述情况下，初始值选为1的时候，从数学上可以证明不存在这两种情况，所以可以放心使用。该方法收敛的速度很快很快，循环十次足够使用，不相信的同志们可以自己循环下看看。）

abel

# Example 3.11 债券价格与必要收益率

票面价值1000，必要收益率从5%到14%，票息率为9%的20年期债券价格

cf <- c(rep(45, 39), 1045)

yeild <- (5:14)/100

(price <- sapply(yeild, function(x) sum(cf*(1/(1+x/2))^(1:40))))

结果为：1502.0555 1346.7216 1213.5507 1098.9639 1000.0000  914.2046  839.5388  774.3055  717.0895  666.7073

plot(yeild, price, type="b", ylim=c(0, 1600))

# Example 3.12 债券价格时间轨迹

面值1000，期限20年，票息率9%，必要收益率12%的债权逼近到期日的债权价格情况。

CouponRate <- 0.09/2

n <- 20*2

yeild <- 0.12/2

(p1 <- (cumsum(1000*CouponRate*(1/(1+yeild))^(1:n)))[seq(2, 40, by=2)])

(p2 <- (1000*(1/(1+yeild))^(1:n))[seq(2, 40, by=2)] )

(p <- p1 + p2)

year <- 1:20

(cbind(year, p1, p2, p))[20:1,]

      year        p1        p2        p

[1,]   20 677.08336  97.22219 774.3055

[2,]   19 668.07086 109.23885 777.3097

[3,]   18 657.94442 122.74077 780.6852

[4,]   17 646.56635 137.91153 784.4779

[5,]   16 633.78195 154.95740 788.7393

[6,]   15 619.41740 174.11013 793.5275

[7,]   14 603.27739 195.63014 798.9075

[8,]   13 585.14248 219.81003 804.9525

[9,]   12 564.76609 246.97855 811.7446

[10,]   11 541.87118 277.50510 819.3763

[11,]   10 516.14645 311.80473 827.9512

[12,]    9 487.24216 350.34379 837.5859

[13,]    8 454.76529 393.64628 848.4116

[14,]    7 418.27428 442.30096 860.5752

[15,]    6 377.27298 496.96936 874.2423

[16,]    5 331.20392 558.39478 889.5987

[17,]    4 279.44072 627.41237 906.8531

[18,]    3 221.27959 704.96054 926.2401

[19,]    2 155.92975 792.09366 948.0234

[20,]    1  82.50267 889.99644 972.4991

plot(year, p, type='b', ylim=c(700, 1000))

当必要收益率为7%的时候：

yeild <- 0.07/2

(p1 <- (cumsum(1000*CouponRate*(1/(1+yeild))^(1:n)))[seq(2, 40, by=2)])

(p2 <- (1000*(1/(1+yeild))^(1:n))[seq(2, 40, by=2)] )

(p <- p1 + p2)

year <- 1:20

(cbind(year, p1, p2, p))[20:1,]

      year        p1       p2        p

[1,]   20 960.97826 252.5725 1213.551

[2,]   19 937.84893 270.5619 1208.411

[3,]   18 913.07222 289.8327 1202.905

[4,]   17 886.53079 310.4761 1197.007

[5,]   16 858.09895 332.5897 1190.689

[6,]   15 827.64204 356.2784 1183.920

[7,]   14 795.01585 381.6543 1176.670

[8,]   13 760.06585 408.8377 1168.904

[9,]   12 722.62654 437.9571 1160.584

[10,]   11 682.52062 469.1506 1151.671

[11,]   10 639.55815 502.5659 1142.124

[12,]    9 593.53568 538.3611 1131.897

[13,]    8 544.23526 576.7059 1120.941

[14,]    7 491.42341 617.7818 1109.205

[15,]    6 434.85005 661.7833 1096.633

[16,]    5 374.24724 708.9188 1083.166

[17,]    4 309.32800 759.4116 1068.740

[18,]    3 239.78489 813.5006 1053.286

[19,]    2 165.28856 871.4422 1036.731

[20,]    1  85.48624 933.5107 1018.997

plot(year, p, type='b')

假定票息率为10，清算日为1997年7月17日，到期日为2003年3月1日，若必要收益率为6.5%，求债券价格轨迹（30/360规则）

# Example 3.14

CouponRate <- 0.10/2

yeild <- 0.065/2

w <- 44/180

time <- w + 0:11

cf <- c(rep(5, 11), 105)

sum(cf * (1+0.0325)^(-time))

例子12和13中出现了必要收益率高于和低于票息率这两种情况，可以看出价格相应的为升高和降低，这算是债券的基本特征之一。

abel

回复5Ldanielbaggio 的问题：

SAS中有GUI方式可以用来做证券计算的，在SAS/ETS模块中也有proc可以来完成，而并非大家看到的一些资料中自己编写代码来完成的。

此外，Matlab中也有Financial toolbox等金融计算相关的工具包若，而且做的相当完善。

我给出的例子，只是用R来一步一步介绍如何计算，未使用R中现成pkgs中的函数，而Rmetrics系列包中，通常需要的功能都有了的。

VaR的back testing在R中早就实现了，你在R资料中搜索下就能知道了——一句话，大部分知道的方法R中已经有现成的了；少部分我们不知道的方法中的绝大部分已经在R中实现了。

abel

<br />
# Data set Files process<br />
path <- "D:/data/SAS_Data/2004/Stoindiv"<br />
fls <- dir(path)<br />
ind <- substr(fls, 1, 2)<br />
shanghai.indv <- fls[ind=="a6"]<br />
shenzhen.indv <- fls[ind=="a0"]<br />
for(fl in shanghai.indv)<br />
    file.copy(from=file.path("D:/data/SAS_Data/2004/Stoindiv", fl), to=file.path("D:/data/SAS_Data/2004/Stoindiv/shanghai", fl))<br />
for(fl in shenzhen.indv)<br />
    file.copy(from=file.path("D:/data/SAS_Data/2004/Stoindiv", fl), to=file.path("D:/data/SAS_Data/2004/Stoindiv/shenzhen", fl))<br />
<br />
# Some SAS procs to transform datasets in SAS into CSV files<br />
<br />
for(fl in shanghai.indv)<br />
    file.remove(file.path("D:/data/SAS_Data/2004/Stoindiv/shanghai", fl))<br />
for(fl in shenzhen.indv)<br />
    file.remove(file.path("D:/data/SAS_Data/2004/Stoindiv/shenzhen", fl))<br />
<br />
<br />
##################################<br />
rtn.log <- function(price){<br />
    x <- price<br />
    r_log <- c(NA, diff(log(x)))<br />
    return(r_log)<br />
}<br />
<br />
ShanghaiPath <- "D:/data/SAS_Data/2004/Stoindiv/shanghai"<br />
ShenzhenPath <- "D:/data/SAS_Data/2004/Stoindiv/shenzhen"<br />
Shanghai.fls <- dir(ShanghaiPath)<br />
Shenzhen.fls <- dir(ShenzhenPath)<br />
<br />
# 上证指数<br />
setwd("D:/data")<br />
a1a0001 <- read.csv("a1a0001.csv", stringsAsFactors=F)[, c("date", "Clpr")]<br />
names(a1a0001) <- c("date", "idx")<br />
a1a0001$date <- as.Date(a1a0001$date)<br />
# 1995年至2000年的数据<br />
ind <- a1a0001$date <= as.Date("2000-12-31") & a1a0001$date >= as.Date("1995-12-29")<br />
a1a0001 <- a1a0001[ind,]    # 共1214个交易日<br />
<br />
setwd(ShanghaiPath)    # 沪市个股处理<br />
sh.rtn <- matrix(NA, nrow=nrow(a1a0001) - 1 , ncol=length(Shanghai.fls))<br />
colnames(sh.rtn) <- Shanghai.fls<br />
rownames(sh.rtn) <- as.character(a1a0001$date[-1])<br />
for(fl in Shanghai.fls){<br />
    dat <- read.csv(fl, stringsAsFactors=F)<br />
    dat$date <- as.Date(dat$date)<br />
    dat <- merge(dat, a1a0001, all.y=T)        # 使用指数来确定交易日，个股会停牌等<br />
    r.log <- rtn.log(dat$Clpr)<br />
<br />
    inds <- apply(dat[, c(4:6, 8:10)], 1, function(x)!any(is.na(x)))<br />
    inds <- (1:nrow(dat))[inds]<br />
    inds <- inds[inds > 1]<br />
<br />
    attach(dat)<br />
    for(ind in inds){<br />
        clpr_adj <- Clpr[ind]*(1 + divrat[ind] + Rigrat[ind] + reisvol[ind]/shrout[ind-1]) <br />
            - rigpr[ind]*Rigrat[ind] - reispr[ind]*reisvol[ind]/shrout[ind-1] + divamt[ind]<br />
        r.log[ind] <- log(clpr_adj) - log(Clpr[ind-1])<br />
    }<br />
    detach()<br />
    sh.rtn[,fl] <- r.log[-1]<br />
}<br />
ind <- apply(sh.rtn, 1, function(x) all(is.na(x)))    # 选出收益率全部是缺失值的日期<br />
sh.rtn <- sh.rtn[!ind, ]<br />
<br />
setwd(ShenzhenPath)    # 深市个股处理<br />
sz.rtn <- matrix(NA, nrow=nrow(a1a0001) - 1 , ncol=length(Shenzhen.fls))<br />
colnames(sz.rtn) <- Shenzhen.fls<br />
rownames(sz.rtn) <- as.character(a1a0001$date[-1])<br />
for(fl in Shenzhen.fls){<br />
    dat <- read.csv(fl, stringsAsFactors=F)<br />
    dat$date <- as.Date(dat$date)<br />
    dat <- merge(dat, a1a0001, all.y=T)        # 使用指数来确定交易日，个股会停牌等<br />
    r.log <- rtn.log(dat$Clpr)<br />
<br />
    inds <- apply(dat[, c(4:6, 8:10)], 1, function(x)!any(is.na(x)))<br />
    inds <- (1:nrow(dat))[inds]<br />
    inds <- inds[inds > 1]<br />
<br />
    attach(dat)<br />
    for(ind in inds){<br />
        clpr_adj <- Clpr[ind]*(1 + divrat[ind] + Rigrat[ind] + reisvol[ind]/shrout[ind-1]) <br />
            - rigpr[ind]*Rigrat[ind] - reispr[ind]*reisvol[ind]/shrout[ind-1] + divamt[ind]<br />
        r.log[ind] <- log(clpr_adj) - log(Clpr[ind-1])<br />
    }<br />
    detach()<br />
    sz.rtn[,fl] <- r.log[-1]<br />
}<br />
ind <- apply(sz.rtn, 1, function(x) all(is.na(x)))    # 选出收益率全部是缺失值的日期<br />
sz.rtn <- sz.rtn[!ind, ]<br />
setwd("D:/data")<br />
<br />
write.csv(sh.rtn, file="ShRtn.csv", row.names=F)<br />
write.csv(sz.rtn, file="SzRtn.csv", row.names=F)<br />
<br />
<br />
<br />
<br />
# 随机选择20只股票，并赋权重，计算组合的收益率<br />
set.seed(12345)<br />
wts <- runif(20, 0, 1)<br />
# wts <- wts/sum(wts) 其实这步不需要<br />
set.seed(12345)<br />
stk.smp <- sh.rtn[, sample(Shanghai.fls, 20)]<br />
port.rtn <- apply(stk.smp, 1, function(x) mean(x, weight=wts, na.rm=T))<br />
ts.plot(port.rtn)<br />
<br />
# 两市合并后再做<br />
rtn <- cbind(sh.rtn, sz.rtn)<br />
# 下面省略若干<br />

从上面的过程不难看出，R批量处理的时候使用循环，远远要比SAS中方便。不过当把所有数据存为SAS的一个datset的时候，可以使用by来一次处理，不过通常遇到的数据文件也许很多个，这样使用R效率就极大提高了。

值得提醒的是，上述代码之所以这种形式，主要是为了和原书思路一致，并未做优化。

BTW：帖子的材料可以随意用，只要注明“转自cos论坛Abel”即可；欢迎大家提出批评意见和建议，尤其欢迎大家讨论金融计算和模型相关的内容，不拘是C++、Matlab、SAS或者其他什么。

从上面的代码不难看出

danielbaggio

多谢！继续对你的帖子予以高度支持！ [s:11]

zuozuo1983

楼主应该写本书《基于R的金融计算建模》

不知道有没有这样的打算

abel

回12L的zuozuo1983，编这样一本入门书倒不用花太多功夫，教学上也用得上。不过开这种课的地方应该不多；书要编的合用，要求还是颇高；先找机会在教学中走一遍看看效果了再考虑出书吧。

pengchy

[quote]一句话，大部分知道的方法R中已经有现成的了；少部分我们不知道的方法中的绝大部分已经在R中实现了。

[/quote]

[s:11][s:11]

abel

<br />
# header.csv与shares.csv<br />
# header:    f0001——股票代码<br />
#        f0002——交易所用代码<br />
#        f0003——股票种类<br />
#        f0004——股票名称<br />
# shares:    f0001——股票代码<br />
#        gb002——股票变动日<br />
#        gb004——总股本<br />
#        A股国家法人股=国家股+境内发起人法人股=gb008+gb009<br />
#        A股总股本=总股本-B股总股本-H股总股本<br />
#        A股流通股=gb016<br />
<br />
header <- read.csv("D:/data/SAS_Data/2004/header_genius.csv")<br />
header <- header[header$F0003=='A' & !is.na(header$F0002), ]<br />
shares <- read.csv("D:/data/SAS_Data/2004/shares_genius.csv")[, c(4, 6, 9, 13, 14, 21, 22, 23)]<br />
<br />
stks <- merge(header, shares)<br />
# Sys.setlocale("LC_TIME", "C")    # 设置LC时间格式，在国内使用的时候，R默认设置并非如此！<br />
<br />
stks$date <- as.Date(substr(stks$Gb002, 1, 9), "%d%b%Y") <br />
<br />
stks.1 <- stks[stks$date <= as.Date("2000-12-31"), ]<br />
stks.1 <- stks.1[order(stks.1$F0001, stks.1$date), ]<br />
<br />
ind <- cumsum(tapply(stks.1$F0001, stks.1$F0001, length))<br />
# ind <- cumsum(as.numeric(table(stks.1$F0001)))<br />
stks.1 <- stks.1[ind, ]<br />
(sums <- sapply(stks.1[, c(6:11)], function(x) sum(as.numeric(x), na.rm=TRUE)))<br />
<br />
(sum.state <- sums["Gb008"] + sums["Gb009"])<br />
(sum.all <- sums["Gb004"] - sums["Gb017"] - sums["Gb018"])<br />
(sum.pub <- sums["Gb016"])<br />
(ratio.state <- sum.state/sum.all)<br />
(ratio.pub <- sum.pub/sum.all)<br />
c(sum.all, sum.state, sum.pub, ratio.state, ratio.pub)<br />
<br />
# 截止2001年8月31日<br />
stks.1 <- stks[stks$date <= as.Date("2001-08-31"), ]<br />
stks.1 <- stks.1[order(stks.1$F0001, stks.1$date), ]<br />
<br />
ind <- cumsum(tapply(stks.1$F0001, stks.1$F0001, length))<br />
# ind <- cumsum(as.numeric(table(stks.1$F0001)))<br />
stks.1 <- stks.1[ind, ]<br />
(sums <- sapply(stks.1[, c(6:11)], function(x) sum(as.numeric(x), na.rm=TRUE)))<br />
<br />
(sum.state <- sums["Gb008"] + sums["Gb009"])<br />
(sum.all <- sums["Gb004"] - sums["Gb017"] - sums["Gb018"])<br />
(sum.pub <- sums["Gb016"])<br />
(ratio.state <- sum.state/sum.all)<br />
(ratio.pub <- sum.pub/sum.all)<br />
c(sum.all, sum.state, sum.pub, ratio.state, ratio.pub)<br />
<br />

# 后续过程省略，因为（1）所需过程和第一个例子雷同，简单重复意义不大；（2）书中所使用代码庞杂，不足取。

# 请注意，在R中，可将上述过程写为函数，进而完成这些过程；通过比较也能发现，R在数据预处理任务上有极大的灵活性和方便性。

abel

<br />
# Chapter5：收益率波动计算<br />
<br />
log.ret <- read.csv("D:/data/SAS_Data/2004/logRet.csv")<br />
<br />
# Simple moving average<br />
sma <- function(x, n=20){<br />
    if(length(x) <= n) return(0)<br />
    rslt <- rep(NA, length(x))<br />
    for(i in 21:length(x)){<br />
        rslt[i] <- sqrt(sum(x[(i-n+1):i])/(n-1))<br />
    }<br />
    return(rslt)<br />
}<br />
<br />
sma(log.ret$rr1a0001)[21:30]<br />
sma(log.ret$rr600652)[21:30]<br />
<br />
# Expoential Weighted Moving Average<br />
ewma <- function(x, w=0.94){<br />
    if(length(x) <= 2) return(0)<br />
    rslt <- rep(NA, length(x))<br />
    rslt[2] <- x[2]<br />
    for(i in 3:length(x)){<br />
        rslt[i] <- w*x[i-1] + (1-w)*x[i]<br />
    }<br />
    return(sqrt(rslt))<br />
}<br />
ewma(log.ret$rr1a0001)[2:10]<br />
ewma(log.ret$rr600652)[2:10]<br />
<br />
# GARCH(1, 1)<br />
require(tseries)<br />
fit1 <- garch(log.ret$r1a0001)<br />
fit2 <- garch(log.ret$r600652)<br />
fit1$fitted.values[, 1]<br />
fit2$fitted.values[, 1]<br />

在R的tseries包非常强悍，建议对时间序列和金融数据处理有兴趣的人重点关注下。

顺便说下对原书的感觉：（1）SAS代码编写思路不够清晰；（2）代码的编写都没indent，及其不方便阅读，到处是省略RUN，而SAS社区建议在一个功能后都加上这个statement，表示其功能的结束；（3）部分内容有错误嫌疑；（4）违背了代码重用原则，导致书中大量功能相同的代码重复出现！

这里建议有兴趣的同志，用其他软件依照原书中的公式计算下，看看是不是有问题！

abel

<br />
# 原书中为ARMA(1, 0)-GARCH(1, 1)型<br />
require(fGarch)<br />
fit1 <- garchFit(~arma(1, 0)+garch(1, 1), data=log.ret$r1a0001)<br />
fit2 <- garchFit(~arma(1, 0)+garch(1, 1), data=log.ret$r600652)<br />
(vol.1 <- fit1@sigma.t)<br />
(vol.2 <- fit2@sigma.t)<br />

这里Abel再给一个建议：

对金融计算和建模有兴趣的人，如果将Rmetrics系列包认真看一遍，并将其中的例子用心揣摩，明白其实质。那么，基本上有九成把握可以说，你可以胜任现在主流的金融计算和建模的工作。

此外，在R的stats包中有一个叫convolve的函数，变化下做可以实现sma，ewma；有兴趣的同志自己试一试，记得我在什么地方写过关于这个函数的一个小介绍的。[s:13]

接着，书中给出了这三种方法计算得到波动率以后作图比较，对R来说，直接使用ts.plot即可，毋庸赘述啦。

书中随后给出了多个股票波动计算的代码，其实在R中，你用一个for循环，一切搞定，所做的不过是把上面的形式写成一个函数而已，如：

for(dat in datas) sma(dat, ...)然后，世界就太平了。

在第三节，介绍了资产组合收益波动率的计算。其实可以分解为两步走：（1）计算组合收益率；（2）计算该收益率的波动率。计算组合收益率的过程在前面章节已经给出了，而这章给出了计算收益率波动率的函数，这样简单的功能相加，就可以得到我们需要的结果了。

abel

<br />
# 多个股票波动率计算; NOTE: 原书中的SAS程序运行需要几十分钟！在R中只需要几分钟。<br />
# （1）使用前面章节计算单只股票日收益率；第二章的程序，保存在SzRtn.csv和ShRtn.csv两个文件中<br />
# （2）计算所有股票的波动率：sma形式，5、10、20日<br />
# （3）计算上述形式下波动率全样本期的均值、最大值、最小值，并输出为数据集<br />
setwd(yourpath)<br />
sh.rtn <- read.csv("ShRtn.csv")<br />
sz.rtn <- read.csv("SzRtn.csv")<br />
rtn <- cbind(sh.rtn, sz.rtn)<br />
colnames(rtn) <- substr(colnames(rtn), 1, 7)<br />
<br />
##########################################################<br />
# (2)<br />
sma.5 <- apply(rtn, 2, <br />
            function(x, width=5) {<br />
                xx <- (x - mean(x, na.rm=T))^2<br />
                sma(xx, n=width) <br />
            })<br />
<br />
sma.10 <- apply(rtn, 2, <br />
            function(x, width=10) {<br />
                xx <- (x - mean(x, na.rm=T))^2<br />
                sma(xx, n=width) <br />
            })<br />
<br />
sma.20 <- apply(rtn, 2, <br />
            function(x, width=20) {<br />
                xx <- (x - mean(x, na.rm=T))^2<br />
                sma(xx, n=width) <br />
            })<br />
write.csv(sma.5, file="SMA5.csv", row.names=F)<br />
write.csv(sma.10, file="SMA10.csv", row.names=F)<br />
write.csv(sma.20, file="SMA20.csv", row.names=F)<br />
<br />
##########################################################<br />
# (3)<br />
sma5.stat <- apply(sma.5, 2, <br />
            function(x) {<br />
                if(!all(is.na(x)))<br />
                    c(mean(x, na.rm=T), max(x, na.rm=T), min(x, na.rm=T))<br />
                else<br />
                    c(NA, NA, NA)<br />
            })<br />
<br />
sma10.stat <- apply(sma.10, 2, <br />
            function(x) {<br />
                if(!all(is.na(x)))<br />
                    c(mean(x, na.rm=T), max(x, na.rm=T), min(x, na.rm=T))<br />
                else<br />
                    c(NA, NA, NA)<br />
            })<br />
<br />
sma20.stat <- apply(sma.20, 2, <br />
            function(x) {<br />
                if(!all(is.na(x)))<br />
                    c(mean(x, na.rm=T), max(x, na.rm=T), min(x, na.rm=T))<br />
                else<br />
                    c(NA, NA, NA)<br />
            })<br />
<br />
#####################################################<br />
# 写为一个函数<br />
fun <- function(dat, width=5){<br />
    xx <- apply(dat, 2, function(x) (x-mean(x, na.rm=TRUE))^2)<br />
    xx <- sma(xx, n=width)<br />
    stat <- apply(xx, 2, function(x){<br />
                        if(all(is.na(x)))<br />
                            c(NA, NA, NA)<br />
                        else<br />
                            c(mean(x, na.rm=TRUE), max(x, na.rm=TRUE), min(x, na.rm=TRUE))<br />
                    })<br />
}<br />
# fun(rtn, 5)<br />
# fun(rtn, 10)<br />
# fun(rtn, 20)<br />
<br />
######################################################################<br />
# 5.3 等权重组合收益波动率<br />
# （1）挑选出1999年前上市的所有A股股票<br />
# （2）计算2000年日百分比收益<br />
# （3）随机抽取一个投资组合<br />
# （4）计算波动率<br />
<br />
ind <- as.Date(rownames(rtn)) < as.Date("2000-01-01")<br />
ind <- apply(rtn[ind, ], 2, function(x) all(is.na(x)))<br />
rtn.2k <- rtn[, !ind]<br />
ind <- (as.Date(rownames(rtn)) > as.Date("2000-01-01")) & (as.Date(rownames(rtn)) < as.Date("2001-01-01"))<br />
rtn.2k <- rtn.2k[ind, ]    # rtn.2k中包含了2000年全年均有交易数据的全部A股股票，并保存这些股票在2000年全年的日收益率<br />
<br />
# 随机抽取20支股票<br />
rtn.smp <- rtn.2k[, sample(1:nrow(rtn.2k), 20)]<br />
rtn.tot <- apply(rtn.smp, 1, function(x) mean(x, na.rm=TRUE))<br />
rtn.xx <- (rtn.tot - mean(rtn.tot))^2<br />
sma(rtn.xx, 5)<br />
sma(rtn.xx, 10)<br />
sma(rtn.xx, 20)<br />
<br />
mat <- matrix(0, ncol=2, nrow=50)<br />
for(i in 1:50){<br />
    rtn.smp <- rtn.2k[, sample(1:nrow(rtn.2k), i*2)]<br />
    rtn.tot <- apply(rtn.smp, 1, function(x) mean(x, na.rm=TRUE))<br />
    mat[i, ] <- c(i*2, sd(rtn.tot))<br />
}<br />
plot(ts(mat[,2], start=1, freq=1/2))<br />

abel

SAS不可谓不强悍，同时也不可谓不麻烦 [s:12]

其实私底下我多次说过朱教授这本SAS书中的代码其实是不够“优雅”的——单单用SAS的语言而言，也有更加优雅的表达方法：

1、SAS宏使用不太规范，有泛滥的嫌疑，且命名随意（不能望名生意）

2、中间结果未予必要的保存——从而导致各个章节之间大量重复代码（尤其其中计算日收益率之类的，几乎每章重复一次）

3、部分代码可优化的空间很大

4、存在若干细微错误，且很难发现——同时我个人感觉SAS中代码中的错误很难发现

5、代码编写（印刷）到书中没有indent，没有每个step结束之后的run等，这些规范应该是良好可读的SAS代码必备的

6、SAS代码编写的时候一大堆——严重缺乏注释，部分缺乏良好的逻辑组织

以上纯属挑刺型的个人意见，故而不推荐大家使用该书学习SAS来处理金融数据；其实SAS在整理好数据以后来建模，倒是我一直喜欢的，因为有丰富的现成的模型以及很多很多的选项和完备的文档。

btw：

1、个人不推荐使用该书学习金融数据处理

2、个人不推荐使用该书学习SAS

3、上述给出章节基本上包括了该书的主要特点，也给出了相应的R代码的解决，余下章节可使用类似方法得到

因此，该书的SAS改编为R代码到此告一段落 [s:13]