名师讲堂｜Stata 最优化与马科维茨有效前沿的实现

本文使用 4 只股票：银泰资源、大悦城、中国天楹和通化金马的历史数据演示了如何使用 Stata 进行股票分析、构建投资组合（夏普比率最大 or 方差最小）及寻找马科维茨有效前沿。

获取数据

附件中提供了四只股票的交易数据：

clear all

*- 使用上证 50 四根票进行组合

use 000975, clear
keep clsprc date
ren clsprc 银泰资源
save a000975, replace

use 000031, clear
keep clsprc date
ren clsprc 大悦城
save a000031, replace

use 000035, clear
keep clsprc date
ren clsprc 中国天楹
save a000035, replace

use 000766, clear
keep clsprc date
ren clsprc 通化金马
save a000766, replace

然后再把这四个文件 merge 到一起：

use a000975, clear
foreach i in "000031" "000035" "000766"{
    merge 1:1 date using a`i'
    drop _m
}
*- 保留 2019 年
keep if inrange(date, mdy(1, 1, 2019), mdy(12, 31, 2019))
*- 或者
*- keep if yofd(date) == 2019
save rawdata, replace

rawdata 数据集是这样的：

list in 1/10

*>     +------------------------------------------------------+
*>     |       date   银泰资源    大悦城       中国天楹     通化金马 |
*>     |------------------------------------------------------|
*>  1. | 2019-01-02      10.07     4.86       5.54       7.01 |
*>  2. | 2019-01-03      10.44     4.86       5.61       6.77 |
*>  3. | 2019-01-04      10.51     4.99       5.56       6.86 |
*>  4. | 2019-01-07      10.36     5.02       5.53       6.92 |
*>  5. | 2019-01-08       9.94        5       5.54       6.97 |
*>     |------------------------------------------------------|
*>  6. | 2019-01-09       9.97     5.04       5.42       6.99 |
*>  7. | 2019-01-10      10.26     4.99       5.39       6.87 |
*>  8. | 2019-01-11       10.2     5.01       5.47       6.88 |
*>  9. | 2019-01-14      10.14        5       5.57        6.8 |
*> 10. | 2019-01-15      10.29     5.06       5.54       6.77 |
*>     +------------------------------------------------------+

观察四只股票的股价走势

把期初的股价标准化为 100，观测四只股票的股价走势关系：

注意：colorscheme 也是外部命令，附件中有其离线安装包。

use rawdata, clear
*- 股价标准化（2019年1月2日的股价设定为100）
foreach i of varlist _all{
    if "`i'" != "date"{
        replace `i' = (`i' / `=`i'[1]') * 100
        format `i' %6.2f
    }
}
colorscheme 4, palette(Dark2)
*- di date("20191231", "YMD")
tw ///
line 银泰资源 大悦城 中国天楹 通化金马 date, ///
    xla(21550(45)21914, ang(20)) ///
    leg(order(1 "银泰资源" 2 "大悦城" ///
        3 "中国天楹" 4 "通化金马") pos(1) ring(0)) ///
    lp(solid solid solid solid) lc(`r(colors)')

接下来计算每只股票的对数收益率，这里需要注意由于股票交易日之间可能存在间隔，所以不能直接用日期作为时间变量：

use rawdata, clear
*- 由于股票交易日之间可能存在间隔，所以不能直接用日期作为时间变量
gen dateid = _n
tsset dateid
foreach i of varlist _all{
    if "`i'" != "date" & "`i'" != "dateid"{
        gen l`i' = l.`i'
        replace `i' = log(`i'/l`i')
        drop l`i'
        format `i' %6.4f
    }
}
drop dateid
save returndata, replace

returndata 是这样的：

list in 1/10

*>     +-------------------------------------------------------+
*>     |       date    银泰资源      大悦城     中国天楹     通化金马 |
*>     |-------------------------------------------------------|
*>  1. | 2019-01-02          .         .          .          . |
*>  2. | 2019-01-03     0.0361    0.0000     0.0126    -0.0348 |
*>  3. | 2019-01-04     0.0067    0.0264    -0.0090     0.0132 |
*>  4. | 2019-01-07    -0.0144    0.0060    -0.0054     0.0087 |
*>  5. | 2019-01-08    -0.0414   -0.0040     0.0018     0.0072 |
*>     |-------------------------------------------------------|
*>  6. | 2019-01-09     0.0030    0.0080    -0.0219     0.0029 |
*>  7. | 2019-01-10     0.0287   -0.0100    -0.0056    -0.0173 |
*>  8. | 2019-01-11    -0.0059    0.0040     0.0147     0.0015 |
*>  9. | 2019-01-14    -0.0059   -0.0020     0.0181    -0.0117 |
*> 10. | 2019-01-15     0.0147    0.0119    -0.0054    -0.0044 |
*>     +-------------------------------------------------------+

观测一下股票对数收益率的分布：

use returndata, clear
foreach i of varlist _all{
    if "`i'" != "date"{
        ren `i' returns`i'
    }
}

reshape long returns, i(date) j(name) string
hist returns, by(name, r note("")) bin(50) ///
    xti(日对数收益率) yti(频率) ///
    xla(, format(%6.2f))  norm ///
    leg(rows(1))
gr export pic2.png, width(1200) replace

这就是“尖峰厚尾”了（峰比正态分布高，尾比正态分布厚）。

投资组合的构建

夏普比率代表交易者每多承担一份风险，相应的就可以拿到几份超额报酬。若其为正，代表回报率高过波动风险；若其为负值，代表操作风险大过回报率。这样，每个投资组合都可以计算夏普比率，即投资回报与多冒风险的比例，这个比率越高对应的投资组合越好。

首先计算每只股票对数收益率的均值，然后存储在矩阵 returns 中（其实是个一维的行向量）：

use returndata, clear
drop in 1

local j = 1
foreach i of varlist _all{
    if "`i'" != "date"{
        qui sum `i'
        local r`j' = r(mean)
        local j = `j' + 1
    }
}
di "`r1', `r2', `r3', `r4'"
*> 0017526058780214, .0016060119467009, .0004432278358542, -.0004437631669057

mat returns = (`r1', `r2', `r3', `r4')

然后计算协方差矩阵，并将协方差矩阵存储在矩阵 variance 中：

corr 银泰资源 大悦城 中国天楹 通化金马, cov

*>         |银泰资源大悦城中国天楹通化金马
*> ---------+------------------------------------
*> 银泰资源  |  .000846
*> 大悦城  |  .000066  .000472
*> 中国天楹  |  .000159  .000181  .000626
*> 通化金马  |  .000165  .000363  .000301  .001137

mat variance = r(C)

接下来就是使用 Stata 进行最优化的，为了便于大家理解，我先举一个 Stata 帮助文档里面的例子。假如我想寻找最大化，代码如下：

mata:
void myeval(todo, x, y, g, H)
{
y = exp(-x^2 + x - 3)
}
S = optimize_init()
optimize_init_evaluator(S, &myeval())
optimize_init_params(S, 0)
x = optimize(S)
x
end

可以看到最优化之后得到的 x 取值是 0.5。我们可以把这个函数的图像画出来看看：

tw function y = exp(-x^2 + x - 3), ///
 yti("y = e{sup:-x{sup:2} + x - 3}") ///
 ylab(, format(%6.2f)) xlab(, format(%6.1f)) ///
 xline(0.5)

gr export pic3.png, width(1200) replace

可以看到，确实是 0.5 处的 y 取最大值。

我再解释一下这段最优化代码，因为它不同于我们常见的 Stata 代码，它是 Mata 的代码：

mata: 表示进入 Mata 语言的环境，然后是一个 myeval 函数，这个函数名前面有个 void，表示这个函数是无返回值的。函数有五个参数 —— todo, x, y, g, H，我们这里只需要用到 x 和 y，剩下三个参数的含义可以自行阅读帮助文档（help optimize）。

也就是说，待调优的参数是放在函数的第二个参数上；函数的返回值是放在第三个参数上。

S = optimize_init() 表示初始化，optimize_init_evaluator(S, &myeval()) 设定要最优化的函数；optimize_init_params(S, 0) 是设定初始值，这里是设定初始的 x = 0，然后默认使用的最优化方法是 Newton–Raphson 方法。最后执行 x = optimize(S) 进行最优化，x 就是我们需要的结果。

我们再回到我们的投资组合构建上。

首先我们允许卖空，就是说权重参数可以为负数：

mata:
mata clear
void sharpe_ratio(real scalar todo, real vector w123, sharpe, g, H)
{
    real scalar w4
    w4 = 1 - sum(w123)
    real vector rmean
    rmean = st_matrix("returns")
    real scalar rp
    real vector w
    w = (w123, w4)
    rp = w * rmean'
    rp = rp * 252
    real matrix variancep
    variancep = st_matrix("variance")
    variancep = w * variancep * w'
    variancep = variancep * 252
    real scalar std
    std = sqrt(variancep[1, 1])
    sharpe = (rp - 0.04) / std
}

S = optimize_init()
optimize_init_evaluator(S, &sharpe_ratio())
optimize_init_technique(S, "nr")
optimize_init_params(S, J(1, 3, 0.25))
wh = optimize(S)
wh

w = (wh, 1-sum(wh))
rp_yearly = (wh, 1-sum(wh)) * st_matrix("returns")' * 252
variancep = (w * st_matrix("variance") * w') * 252
std = sqrt(variancep[1, 1])
sharpe = (rp_yearly - 0.04) / std
"权重分别为："
w

*>                 1              2              3              4
*>   +-------------------------------------------------------------+
*> 1 |   .4940084988    1.125948699   -.0810411934   -.5389160044  |
*>   +-------------------------------------------------------------+

"此时的组合收益为："
rp_yearly

*> .7250847964

"此时的组合标准差为："
std

*> .4095752946

"此时组合的夏普比率为："
sharpe

*>  1.672671192

end

这就是允许卖空时夏普比率最大的组合。

再考虑不允许卖空的情形，为了保证权重为正且总和为 1，我借用了函数的特点：

use returndata, clear
drop in 1

*- 计算每只股票对数收益率的均值
local j = 1
foreach i of varlist _all{
    if "`i'" != "date"{
        qui sum `i'
        local r`j' = r(mean)
        local j = `j' + 1
    }
}
di "`r1', `r2', `r3', `r4'"
mat returns = (`r1', `r2', `r3', `r4')

*- 计算对数收益率的协方差矩阵
corr 银泰资源 大悦城 中国天楹 通化金马, cov
mat variance = r(C)

mata:
mata clear
void sharpe_ratio_no_short(real scalar todo, real vector w1234, sharpe, g, H)
{

    real vector rmean
    rmean = st_matrix("returns")
    real scalar rp
    real vector w
    w = (exp(w1234[1]), exp(w1234[2]), exp(w1234[3]), exp(w1234[4])) / sum((exp(w1234[1]), exp(w1234[2]), exp(w1234[3]), exp(w1234[4])))
    rp = w * rmean'
    rp = rp * 252
    real matrix variancep
    variancep = st_matrix("variance")
    variancep = w * variancep * w'
    variancep = variancep * 252
    real scalar std
    std = sqrt(variancep[1, 1])
    sharpe = (rp - 0.04) / std
}

S = optimize_init()
optimize_init_evaluator(S, &sharpe_ratio_no_short())
optimize_init_technique(S, "nr")
optimize_init_params(S, J(1, 4, 0.01))
wh = optimize(S)
wh

w = (exp(wh[1]), exp(wh[2]), exp(wh[3]), exp(wh[4])) / sum((exp(wh[1]), exp(wh[2]), exp(wh[3]), exp(wh[4])))
rp_yearly = ((exp(wh[1]), exp(wh[2]), exp(wh[3]), exp(wh[4])) / sum((exp(wh[1]), exp(wh[2]), exp(wh[3]), exp(wh[4])))) * st_matrix("returns")' * 252
variancep = (w * st_matrix("variance") * w') * 252
std = sqrt(variancep[1, 1])
sharpe = (rp_yearly - 0.04) / std
"权重分别为："
w

*>                 1             2             3             4
*>    +---------------------------------------------------------+
*>  1 |  .3700694896   .6299305104   6.33387e-22   3.50516e-18  |
*>    +---------------------------------------------------------+

"此时的组合收益为："
rp_yearly

*> .4183859958

"此时的组合标准差为："
std

*> .2901268944

"此时组合的夏普比率为："
sharpe

*> 1.304208617

end

不过需要注意的是，这里的最优化结果并没有达到收敛。能不能收敛跟数据是有关系的，无法收敛就说明在不允许卖空的时候无法达到夏普比率的最优化。注意到权重，实际上这个结果是让我们只买前两只股票。只要买一点点后两只股票，夏普比率就会变小。

夏普比率函数

我还写了个可以计算夏普比率的函数：

*- 夏普比率函数
use returndata, clear
drop in 1

*- 计算每只股票对数收益率的均值
local j = 1
foreach i of varlist _all{
    if "`i'" != "date"{
        qui sum `i'
        local r`j' = r(mean)
        local j = `j' + 1
    }
}
di "`r1', `r2', `r3', `r4'"
mat returns = (`r1', `r2', `r3', `r4')

*- 计算对数收益率的协方差矩阵
corr 银泰资源 大悦城 中国天楹 通化金马, cov
mat variance = r(C)
mata:
mata clear
real scalar sharpe_ratio(real scalar w1, real scalar w2, real scalar w3)
{
    real scalar w4
    w4 = 1 - w1 - w2 - w3
    real colvector w
    w = (w1 \ w2 \ w3 \ w4)
    real matrix r
    real colvector rmean
    rmean = st_matrix("returns")
    real scalar rp
    rp = w' * rmean'
    rp = rp * 252
    real matrix variancep
    variancep = st_matrix("variance")
    variancep = w' * variancep * w
    variancep = variancep * 252
    real scalar std
    std = sqrt(variancep[1, 1])
    real scalar sharpe
    sharpe = (rp - 0.04) / std
    return(sharpe)
}
end

mata: sharpe_ratio(0.5, 0.5, 0)

结果：

mata: sharpe_ratio(0.3700694896, 1-0.3700694896, 0)
*> 1.304208617

和上面最优化运行得到的结果是一致的，也就是说，在不允许卖空的时候，夏普比率最大的组合是买入 0.37 份的银泰资源和 0.63 份的大悦城，剩下两只股票不买入。

下面我们再用蒙特卡洛模拟观测不允许卖空时，所有可能的组合的收益和标准差的关系：

cap prog drop front
prog def front, rclass
    version 15.1
    use returndata, clear
    local w1 = runiform()
    local w2 = runiform()
    local w3 = runiform()
    local w4 = runiform()
    mat weight = (`w1' \ `w2' \ `w3' \ `w4')
    mat list weight
    mat weight = weight / (`w1'+`w2'+`w3'+`w4')
    mat list weight
    ret scalar w1 = weight[1, 1]
    ret scalar w2 = weight[2, 1]
    ret scalar w3 = weight[3, 1]
    ret scalar w4 = weight[4, 1]
    local j = 1
    foreach i of varlist _all{
        if "`i'" != "date"{
            qui sum `i'
            local r`j' = r(mean)
            local j = `j' + 1
        }
    }
    mat returns = (`r1', `r2', `r3', `r4')
    mat list returns
    mat a =  returns * weight * 252
    mat list a
    ret scalar ret = a[1, 1]
    corr 银泰资源 大悦城 中国天楹 通化金马, cov
    ret list
    mat variance = r(C)
    mat list variance
    mat b =  weight' * variance * 252 * weight
    mat list b
    ret scalar var = b[1, 1]
    ret scalar std = sqrt(b[1, 1])
end

simulate ret = r(ret) var = r(var) std = r(std) w1 = r(w1) w2 = r(w2) w3 = r(w3) w4 = r(w4), reps(100000): front
*- 计算夏普比率
gen sharpe_ratio = (ret - 0.04)/std
save montecarlo, replace

绘图展示：

use montecarlo, clear
tw ///
sc ret std, msize(*0.01) xti(波动率) yti(收益率) msymbol(o) ///
    xla(#6, format(%6.2f)) yla(, format(%6.1f)) ///
    ti("投资组合收益率与波动率的关系", size(*1.2)) || ///
scatteri 0.7250847964 0.4095752946 "允许卖空", mc("red") msymbol(o) || ///
scatteri 0.4183859958 0.2901268944 "不允许卖空", mc("green") msymbol(o) ///
    mlabp(2) leg(off)

gr export pic4.png, width(1200) replace

显然，卖空即提高了收益又增加了风险。

最小化方差

最小化方差也就是最小化标准差，显然就是上图中最右边的那个点了。下面用程序寻找这个点对于的组合构建方法。

首先还是允许做空：

use returndata, clear
drop in 1
local j = 1
foreach i of varlist _all{
    if "`i'" != "date"{
        qui sum `i'
        local r`j' = r(mean)
        local j = `j' + 1
    }
}
di "`r1', `r2', `r3', `r4'"
mat returns = (`r1', `r2', `r3', `r4')
corr 银泰资源 大悦城 中国天楹 通化金马, cov
mat variance = r(C)
mata:
mata clear
void min_std(real scalar todo, real vector w123, std, g, H)
{
    real scalar w4
    w4 = 1 - sum(w123)
    real vector rmean
    rmean = st_matrix("returns")
    real scalar rp
    real vector w
    w = (w123, w4)
    rp = w * rmean'
    rp = rp * 252
    real matrix variancep
    variancep = st_matrix("variance")
    variancep = w * variancep * w'
    variancep = variancep * 252
    std = -sqrt(variancep[1, 1])
}

S = optimize_init()
optimize_init_evaluator(S, &min_std())
optimize_init_technique(S, "nr")
optimize_init_params(S, J(1, 3, 0.25))
wh = optimize(S)
wh

w = (wh, 1-sum(wh))
rp_yearly = (wh, 1-sum(wh)) * st_matrix("returns")' * 252
variancep = (w * st_matrix("variance") * w') * 252
std = sqrt(variancep[1, 1])
sharpe = (rp_yearly - 0.04) / std
"权重分别为："
w

"此时的组合收益为："
rp_yearly

*> 0.3399958417

"此时的组合标准差为："
std

*> 0.271031146

"此时组合的夏普比率为："
sharpe

*> 1.106868513

end

注意上面的代码里面，我在 std = -sqrt(variancep[1, 1]) 里面加了个负号，这样最大化 std 的时候就最小化标准差了。

这个点在这里：

use montecarlo, clear
tw ///
sc ret std, msize(*0.01) xti(波动率) yti(收益率) msymbol(o) ///
    xla(#6, format(%6.2f)) yla(, format(%6.1f)) ///
    ti("投资组合收益率与波动率的关系", size(*1.2)) || ///
scatteri 0.7250847964 0.4095752946 "允许卖空", ///
 mc("red") msymbol(o) || ///
scatteri 0.4183859958 0.2901268944 "不允许卖空", ///
 mc("green") msymbol(o) ///
    mlabp(2) leg(off) || ///
scatteri 0.3399958417 0.271031146 "最小方差点", ///
    mc("pink")  msymbol(o) mlabc("pink")
gr export pic5.png, width(1200) replace

构造有效前沿

有效前沿是一条线，其上的点表示既定收益率下方差最小的投资组合。显然就是最小方差点到最大夏普比率点的边界线。从图中可以看出从最小方差到允许卖空的最大夏普比率点的收益率范围为：[0.33, 0.73]，以 0.01 为间隔，计算每个收益率水平上的最小方差。

use returndata, clear
drop in 1
local j = 1
foreach i of varlist _all{
    if "`i'" != "date"{
        qui sum `i'
        local r`j' = r(mean)
        local j = `j' + 1
    }
}
di "`r1', `r2', `r3', `r4'"
mat returns = (`r1', `r2', `r3', `r4')
corr 银泰资源 大悦城 中国天楹 通化金马, cov
mat variance = r(C)
*- 允许卖空
mata:
mata clear
void min_std(real scalar todo, real vector w1234, std, g, H)
{
    real vector rmean
    rmean = st_matrix("returns")
    real scalar rp
    real vector w
    w = w1234
    rp = w * rmean'
    rp = rp * 252
    real matrix variancep
    variancep = st_matrix("variance")
    variancep = w * variancep * w'
    variancep = variancep * 252
    std = -sqrt(variancep[1, 1])
}
std_vector = J(1, 41, 0)
j = 1
for(i = 0.33; i < 0.73; i = i + 0.01)
{
    S = optimize_init()
    optimize_init_evaluator(S, &min_std())
    optimize_init_technique(S, "nr")
    optimize_init_constraints(S, ((st_matrix("returns") * 252 \ J(1, 4, 1)), (i \ 1)))
    optimize_init_params(S, J(1, 4, 0.25))
    wh = optimize(S)
    wh

    w = wh
    variancep = (w * st_matrix("variance") * w') * 252
    std = sqrt(variancep[1, 1])
    std_vector[1, j] = std
    j++
}
st_matrix("std_vector", std_vector)
end
clear
set obs 41
gen ret = (_n + 32) / 100
gen std = .
forval i = 1/41{
    replace std = std_vector[1, `i'] in `i'
}
gen front = 1
save front, replace

这里我使用的是含约束的最优化，通过设定 optimize_init_constraints() 的约束矩阵实现的。

同样我们把刚刚得到的 41 个点（得删除 0.73 对应的点，这个点的 std = 0，因为允许卖空时的最大收益才 0.725）绘图展示：

use montecarlo, clear
append using front
drop if std == 0

tw ///
sc ret std, msize(*0.01) xti(波动率) yti(收益率) msymbol(o) ///
    xla(#6, format(%6.2f)) yla(, format(%6.1f)) ///
    ti("投资组合收益率与波动率的关系", size(*1.2)) || ///
scatteri 0.7250847964 0.4095752946 "允许卖空", ///
 mc("red") msymbol(o) || ///
scatteri 0.4183859958 0.2901268944 "不允许卖空", ///
 mc("green") msymbol(o) ///
    mlabp(2) leg(off) || ///
scatteri 0.3399958417 0.271031146 "最小方差点", ///
    mc("pink")  msymbol(o) mlabc("pink") || ///
sc ret std if front == 1, msize(tiny) msymbol(D) ///
    mc("dkorange") leg(off)
gr export pic6.png, width(1200) replace

感觉用 Stata 做金融分析也是蛮有意思的。

点击这里跳转到 RStata 短书平台获取附件：名师讲堂｜Stata 最优化与马科维茨有效前沿的实现