2026a

# 基于求解器的投资组合优化问题的二次规划


说明关于基本投资组合模型的基于求解器的二次规划的示例。

此示例说明如何使用 quadprog 中的内点二次规划算法来求解投资组合优化问题。函数 quadprog 属于 Optimization Toolbox。

在此示例中定义问题的矩阵是稠密矩阵;然而,quadprog 中的内点算法也可以利用问题矩阵的稀疏性来提高速度。

# 二次模型

假设存在 种不同资产。资产回报率 是随机变量,预期值为 。问题是在指定的最低预期回报率约束下,找出在每项资产 上的投资比例 为多少,才能将风险降至最低。

表示资产回报率的协方差矩阵。

经典的均值-方差模型是在一组约束条件下最小化投资组合风险,风险可由下式衡量:

其中约束如下:

  • 预期回报率不应低于投资者期望的最低投资组合回报率

  • 投资组合比例 的总和应为 1,

  • 而且作为比例(或百分比),它们应为介于 0 和 1 之间的数字,

由于最小化投资组合风险的目标函数是二次型的,而约束是线性的,因此该优化问题是二次规划,即 QP。

# 225 项资产问题

现在让我们求解包含 225 项资产的 QP。该数据集来自 OR-Library [Chang, T.-J., Meade, N., Beasley, J.E. and Sharaiha, Y.M., "Heuristics for cardinality constrained portfolio optimisation" Computers & Operations Research 27 (2000) 1271-1302]。

我们加载数据集,然后按 quadprog 所需的格式设置约束。在此数据集中,回报率 的范围在 -0.008489 和 0.003971 之间;我们在两者之间选择一个期望的回报率 ,例如 0.002 (0.2%)。

加载存储在 MAT 文件中的数据集。

using TyOptimization
using TyMath
using TyPlot
using TyBase
path = pkgdir(TyOptimization)
data = load(path * "/data/port5.mat")

基于相关矩阵计算协方差矩阵。

Covariance = Correlation .* (stdDev_return * stdDev_return');
nAssets = length(mean_return); r = 0.002;     # number of assets and desired return
Aeq = ones(1,nAssets); beq = 1;              # equality Aeq*x = beq
Aineq = -mean_return'; bineq = -r;           # inequality Aineq*x <= bineq
lb = zeros(nAssets,1); ub = ones(nAssets,1); # bounds lb <= x <= ub
c = zeros(nAssets,1);                        # objective has no linear term; set it to zero

# 在 Quadprog 中选择内点算法

为了使用内点算法求解 QP,我们将选项算法设置为 'interior-point-convex'。

options = optimoptions(:quadprog,Algorithm = "interior-point-convex");

# 求解 225 项资产问题

我们现在设置一些额外的选项,并调用求解器 quadprog。

设置附加选项:打开迭代输出,并设置更严格的最优性终止容差。

options = optimoptions(options,Display = "iter");

调用求解器并测量挂钟时间。

@time x1,fval1 = quadprog(Covariance,vec(c),Aineq,[bineq],Aeq,[beq],vec(lb),vec(ub),[],options);
Iter            Fval   Primal Infeas     Dual Infeas   Complementarity
   0    7.212813e+00    1.975147e+03    2.230265e-01      2.121324e+01
   1    5.087187e-04    1.360917e+00    1.536699e-04      2.830551e+00
   2    3.936395e-04    8.181869e-01    9.238672e-05      4.511132e+00
   3    3.525742e-04    2.468966e-01    2.787868e-05      7.635321e-01
   4    5.519330e-04    1.234483e-05    1.393936e-09      1.487781e+00
   5    5.519320e-04    1.923520e-09    2.455588e-13      2.376074e-04
   6    4.680051e-04    3.926649e-10    5.012570e-14      7.796372e-05
   7    3.606961e-04    1.377623e-11    1.753372e-15      1.406325e-05
   8    2.982304e-04    1.989739e-14    3.686287e-18      1.085233e-05
   9    2.729816e-04    1.874721e-14    2.493665e-18      4.602483e-06
  10    2.238534e-04    1.014434e-14    5.854692e-17      4.445010e-05
  11    2.099378e-04    1.628271e-14    2.818926e-18      6.357395e-06
  12    1.971963e-04    1.986036e-14    1.246832e-18      7.627446e-07
  13    1.952795e-04    1.820585e-14    8.673617e-19      3.467287e-08
  14    1.949127e-04    1.989459e-14    1.002887e-18      1.305679e-10

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance and constraints are satisfied to within the value of the constraint tolerance.

  0.282460 seconds (5.12 k allocations: 64.743 MiB, 12.99% gc time)

绘制结果。

bar(x1,width = 2)
ylabel("Fraction of investment")
xlabel("Assets")
title("Standard model - 225-asset problem")
grid("on")

# 具有组约束的 225 项资产问题

我们现在向模型增加一组约束,即要求投资者 30% 的资金必须投资于资产 1 至 75,30% 投资于资产 76 至 150,30% 投资于资产 151 至 225。每组资产可以是不同行业,例如技术、汽车和制药。

实现此新要求的约束为:

向现有等式约束添加这组不等式约束。

Groups = blkdiag(ones(1,div(nAssets,3)),ones(1,div(nAssets,3)),ones(1,div(nAssets,3)));
Aineq = [Aineq; -Groups];         # convert to <= constraint
bineq = [bineq; -0.3*ones(3,1)];  # by changing signs

调用求解器并测量挂钟时间。

@time x2,fval2 = quadprog(Covariance,vec(c),Aineq,vec(bineq),Aeq,[beq],vec(lb),vec(ub),[],options);
Iter            Fval   Primal Infeas     Dual Infeas   Complementarity
   0    7.212813e+00    2.304305e+03    5.062751e+00      2.598079e+01
   1    5.014626e-04    1.446136e+00    3.177280e-03      3.299995e+00
   2    4.043363e-04    8.119089e-01    1.783832e-03      2.442291e+00
   3    3.474328e-04    2.375218e-01    5.218552e-04      1.628265e+00
   4    3.829175e-04    1.187609e-05    2.609276e-08      3.536073e-01
   5    3.829158e-04    1.191021e-09    2.628376e-12      3.943413e-05
   6    3.536513e-04    2.790679e-10    6.158524e-13      1.232625e-05
   7    2.965677e-04    2.469723e-11    5.447932e-14      5.401752e-06
   8    2.621091e-04    1.974690e-14    7.209944e-18      2.979064e-06
   9    2.234785e-04    1.310451e-14    2.276825e-18      5.842570e-06
  10    2.082052e-04    1.542020e-14    1.084202e-18      1.104407e-06
  11    2.027746e-04    1.886318e-14    2.114194e-18      3.336141e-08
  12    2.009830e-04    2.021292e-14    9.486769e-19      1.646531e-09

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance and constraints are satisfied to within the value of the constraint tolerance.

  0.240488 seconds (5.03 k allocations: 62.282 MiB, 6.91% gc time)

绘制结果,叠加在先前问题的结果上。

hold("on")
bar(x2)
legend(["Standard","With group constranint"])
hold("off")

# 迄今为止结果摘要

我们从第二个条形图中看到,由于存在额外的组约束,与第一个投资组合相比,现在投资组合在三个资产组中的分布更均匀。这种强加的多样化也导致风险略微增加,这是通过目标函数来测量的(请参阅两次运行的迭代输出中最后一次迭代的标记为“f(x)”的列)。

# 使用随机数据的 1000 项资产问题

为了说明 quadprog 的内点算法如何处理更大型的问题,我们将使用包含 1000 项资产的随机生成的数据集。

重置随机流以实现可再现性。

rng = MT19937ar(0);
nAssets = 1000; # desired number of assets

生成 -0.1 到 0.4 之间的回报率均值数据。

a = -0.1; b = 0.4;
mean_return = a .+ (b-a).*rand(rng,nAssets,1);

生成 0.08 到 0.6 之间的回报率标准差数据。

a = 0.08; b = 0.6;
stdDev_return = a .+ (b-a).*rand(rng,nAssets,1);
# Correlation matrix, generated using Correlation = gallery('randcorr',nAssets).
# (Generating a correlation matrix of this size takes a while, so we load
# a pre-generated one instead.)
path = pkgdir(TyOptimization)
data = load(path * "/data/correlationMatrixDemo.mat")
# Calculate covariance matrix from correlation matrix.
Covariance = Correlation .* (stdDev_return * stdDev_return');

# 定义并求解随机生成的 1000 项资产问题

我们现在定义标准的 QP 问题(此处没有组约束)并求解。

r = 0.15;                                     # desired return
Aeq = ones(1,nAssets); beq = [1];               # equality Aeq*x = beq
Aineq = -mean_return'; bineq = [-r];            # inequality Aineq*x <= bineq
lb = zeros(nAssets); ub = ones(nAssets);  # bounds lb <= x <= ub
c = zeros(nAssets);                         # objective has no linear term; set it to zero

@time x3, = quadprog(Covariance,c,Aineq,bineq,Aeq,beq,lb,ub,[],options);
Iter            Fval   Primal Infeas     Dual Infeas   Complementarity
   0    2.144124e+01    1.306964e+04    2.686312e+00      4.519961e+01
   1    9.603334e-02    1.156196e+02    2.376427e-02      5.507976e+00
   2    9.472022e-05    9.980635e-01    2.051404e-04      9.096405e-02
   3    7.912080e-05    4.990318e-05    1.025702e-08      1.035340e-03
   4    7.905810e-05    2.067984e-08    4.250529e-12      1.123253e-06
   5    3.353342e-05    2.173559e-10    4.471614e-14      4.747129e-07
   6    1.856073e-05    6.550039e-14    1.111307e-17      4.743841e-07
   7    1.458766e-05    8.263706e-14    2.289530e-18      2.539410e-08
   8    1.371413e-05    8.318969e-14    9.911344e-19      9.405355e-10

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance and constraints are satisfied to within the value of the constraint tolerance.

  2.605711 seconds (15.60 k allocations: 722.423 MiB, 2.48% gc time)