不借助 matlab 内置函数，生撸均值方差模型

前言

我在之前的一篇文章中介绍了，如何使用 matlab 自带的函数（对象）portfolio,实现均值方差模型。matlab 内置的函数自然实用。但是对于均值方差这样的基本模型，我们不能总是依赖 matlab ，调包侠是没前途的，要努力成为造轮子的人。

具体要求

这其实是一个课程作业，接下来我们看看这个作业有什么要求？

当然这个作业比较有趣啊，我们不仅不能用 portfolio 对象，连 cov,std 之类的函数都不能使用，得自己编写。至于本文需要用到的数据和代码，我在文末附有连接。

均值方差模型

这个作业中，算标准差、协方差、几何收益率之类都是简单的小问题，比较难以实现的是均值方差模型。

均值方差模型就是要求解这个优化问题，图中的向量 x 是投资组合权重。Q 是协方差，A 是资产的期望收益率，b 是投资组合期望收益率。x 在0到1之类意味着这里不考虑卖空操作，如果有卖空，把 lower bounds 改为负无穷即可。

而要解这个问题，我们实际上还是不得不求助于 matlab Optimization Toolbox 工具箱中的 solver。其实但就这个问题而言，只是用到了 quadprog 函数，该函数的详细用法参考 matlab 帮助文档即可。我下面的代码中其实也给了比较详细的注释，耐心点应该完全可以看懂。

手撸代码

接下来，就是激动人心的部分。

计算月收益率

首先，我们需要编写函数计算月收益率。注意这个月收益率，老师要求是这样算的：

To compute a historical monthly return, use the closing
price for the first trading day of a month and the closing price of the last trading day of the month.

function ret = calc_mon_ret(f)
% get month return

% 从 csv 文件中读取数据
[num,txt] = xlsread(f);

adj_close = num(:,5);
trade_date = txt(2:end,1);

date_vec = datevec(trade_date);
% 把日数据的年和月提取出来，这么写其实麻烦了，不过也懒得改了
year = unique(date_vec(:,1));
month = unique(date_vec(:,2));

ret = [];
for i = 1:length(year)
    for j = 1:length(month)
    % temp 就是同一个月的数据
        temp = date_vec(:,1) == year(i) & date_vec(:,2) == month(j);
        if any(temp)
            day = find(temp);
            % 计算收益率
            temp_ret = adj_close(day(end))/adj_close(day(1)) - 1;
            ret = [ret;temp_ret];
        end
    end
end

end

计算协方差

其实小问题中就协方差计算有点难度。

function value = my_cov(retmat)
% caculate covariation

[m,n] = size(retmat);
value = nan(n,n);
% arithmetic mean
am = sum(retmat) / m;

for i = 1:n
    for j = 1:n
        temp_i = retmat(:,i) - am(i);
        temp_j = retmat(:,j) - am(j);
        temp = sum(temp_i .* temp_j);
        value(i,j) = temp / m;
    end
end

end

函数方面，我觉得自己写的没有问题，但是这个结果和 matlab 自带函数 cov 计算出来有些许差别，我不知道原因在哪，有知道的朋友可以指正一下。

均值方差模型的计算

稍微复杂的函数，已经被我们分离成子函数了，接下来该写这个作业的主要脚本了。

close all
clear
clc
%% load data
f1 = 'SPY(Daily).csv';
f2 = 'GOVT(Daily).csv';
f3 = 'EEMV(Daily).csv';

f4 = 'CME(Daily).csv';
f5 = 'BR(Daily).csv';
f6 = 'CBOE(Daily).csv';
f7 = 'ICE(Daily).csv';
f8 = 'ACN(Daily).csv';

variable = {'SPY','GOVT','EEMV','CME','BR','CBOE','ICE','ACN'};

% caculate monthly return
SPY_ret = calc_mon_ret(f1);
GOVT_ret = calc_mon_ret(f2);
EEMV_ret = calc_mon_ret(f3);

CME_ret = calc_mon_ret(f4);
BR_ret = calc_mon_ret(f5);
CBOE_ret = calc_mon_ret(f6);
ICE_ret = calc_mon_ret(f7);
ACN_ret = calc_mon_ret(f8);
% ret matrix
retmat = [SPY_ret GOVT_ret EEMV_ret CME_ret BR_ret CBOE_ret ICE_ret ACN_ret];
%%  mean
% arithmetic mean
[m,n] = size(retmat);
am = sum(retmat) / m;
% Geometric mean
gm = cumprod(retmat + 1);
gm = gm(end,:);
% how many month
nMon = size(retmat,1);
gm = power(gm,1/nMon) - 1;

%% standard deviation
% sd = std(retmat);
sd = nan(1,n);
for i = 1:n
    temp = retmat(:,i) - am(i);
    sd(i) = sqrt( sum(temp .^ 2) / (m-1));
end
%% co-variance
co_variance = my_cov(retmat);
% write co-variance 
temp = num2cell(co_variance);
temp = [variable;temp];
para = [{'covariance'};variable'];
temp = [para temp];
xlswrite('covariance.xls',temp);
%% write parameter to xls
value = [am;gm;sd];
value = num2cell(value);
para = {'';'Arithmetic mean';'Geometric mean';'Standard deviation'};
temp = [variable;value];
temp = [para temp];
xlswrite('mean and sd.xls',temp);
%%  part a MVO

% part a portfolio
retmat_a = retmat(:,1:3);
% parameter
Q = my_cov(retmat_a);
f = zeros(1,3);
A = gm(1:3);
% expected return
b = 0.0001:0.0004:0.008;
b = b';
Aeq = ones(1,3);
beq = 1;
% without short sale
lb = zeros(1,3);
ub = ones(1,3);

% n iteration
n = length(b);
 %weights distribution
w_a = nan(n,3);
%volatility
risk_a = nan(n,1);
 
%iteration
for i = 1:n
    [w_a(i,:),var]=quadprog(Q,f,-A,-b(i),Aeq,beq,lb,ub);
    % volatility
    risk_a(i)= sqrt(var); 
 
end
 
% plot
p = plot(risk_a,b,'g-o');
set(p,'LineWidth',1.5,'MarkerSize',3)
title('The efficient frontier of MVO(three assets)');
xlabel('volatility');
ylabel('Expected Monthly Return');
saveas(gcf,'part_a frontier.jpg');
% save data
v = [variable(1:3) {'Expect return','Volatility'}];
value = [w_a b risk_a];
value = num2cell(value);
temp = [v;value];
xlswrite('part a weight and expect goal.xls',temp);

%% part b MVO

% parameter
nAsset = 8;
Q = my_cov(retmat);
f = zeros(1,nAsset);
A = gm;
Aeq = ones(1,nAsset);
beq = 1;
% without short sale
lb = zeros(1,nAsset);
ub = ones(1,nAsset);

 %weights distribution
w_b = nan(n,nAsset);
%volatility
risk_b = nan(n,1);
 
%iteration
for i = 1:n
    [w_b(i,:),var]=quadprog(Q,f,-A,-b(i),Aeq,beq,lb,ub);
    % volatility
    risk_b(i)= sqrt(var); 
end
 
% plot
hold on
p1 = plot(risk_b,b,'r-o');
set(p1,'LineWidth',1.5,'MarkerSize',3)
title('The efficient frontier of MVO');
xlabel('volatility');
ylabel('Expected Monthly Return');
legend('three assets','eight assets','Location','best');
saveas(gcf,'part_b frontier.jpg');

% save data
v = [variable {'Expect return','Volatility'}];
value = [w_b b risk_b];
value = num2cell(value);
temp = [v;value];
xlswrite('part b weight and expect goal.xls',temp);

这部分要提的就是 quadprog 函数的运用，我们通过循环求解出了每一个期望投资组合收益率下的投资组合权重还有标准差。

代码和资料下载

本文代码和相关资料已经上传至 **** 资源下载页面，欢迎下载。如有错误，请不吝赐教。