我在之前的一篇文章中介绍了，如何使用 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 value = my_cov(retmat)% caculate covariation[m,n] = size(retmat);value = nan(n,n);% arithmetic meanam = 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; endendend函数方面，我觉得自己写的没有问题，但是这个结果和 matlab 自带函数 cov 计算出来有些许差别，我不知道原因在哪，有知道的朋友可以指正一下。
均值方差模型的计算稍微复杂的函数，已经被我们分离xxdzxc函数了，接下来该写这个作业的主要脚本了。
close allclearclc%% load dataf1 = '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 returnSPY_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 matrixretmat = [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 meangm = cumprod(retmat + 1);gm = gm(end,:);% how many monthnMon = 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-varianceco_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 xlsvalue = [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 portfolioretmat_a = retmat(:,1:3);% parameterQ = my_cov(retmat_a);f = zeros(1,3);A = gm(1:3);% expected returnb = 0.0001:0.0004:0.008;b = b';Aeq = ones(1,3);beq = 1;% without short salelb = zeros(1,3);ub = ones(1,3);% n iterationn = length(b); %weights distributionw_a = nan(n,3);%volatilityrisk_a = nan(n,1); %iterationfor i = 1:n [w_a(i,:),var]=quadprog(Q,f,-A,-b(i),Aeq,beq,lb,ub); % volatility risk_a(i)= sqrt(var); end % plotp = 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 datav = [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% parameternAsset = 8;Q = my_cov(retmat);f = zeros(1,nAsset);A = gm;Aeq = ones(1,nAsset);beq = 1;% without short salelb = zeros(1,nAsset);ub = ones(1,nAsset); %weights distributionw_b = nan(n,nAsset);%volatilityrisk_b = nan(n,1); %iterationfor i = 1:n [w_b(i,:),var]=quadprog(Q,f,-A,-b(i),Aeq,beq,lb,ub); % volatility risk_b(i)= sqrt(var); end % plothold onp1 = 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 datav = [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 函数的运用，我们通过循环求解出了每一个期望投资组合收益率下的投资组合权重还有标准差。
代码和资料下载本文代码和相关资料已经上传至 csdn 资源下载页面，欢迎下载。如有错误，请不吝赐教。