内容是使用Jupyter Notebook写的，它比较便利编写公式和调试代码．这里的内容是粘贴的，所以代码的排版可能存在无法直接执行．
金融中的线性问题
传统的投资组合管理中，资产配置通常遵循线性模型。
我们可以将投资组合配置问题转换为包含等式和不等式的线性方程组：
$$Ax=B$$
其中，A为已知系数，B为观测值，$x$为未知向量。$x$代表最优的资产配比，可以利用线性代数中的直接或间接方法快速求解。
2.1资本资产定价模型与证券市场线 capital asset pricing model (CAPM)
证券资产风险与回报率关系：
$$R_i = R_f + \beta_i (R_{mkt} - R_f)$$
对于证券$i$，设其回报率为$R_i$，$\beta$系数为$\beta_i$，则回报率$R_i$等于无风险回报率$R_f$与$\beta$系数市场风险溢价之和。
from scipy import stats
stock_returns = [0.065, 0.0265, -0.0593, -0.001, 0.0346]
mkt_returns = [0.055, -0.09, -0.041, 0.045, 0.022]
beta, alpha, r_value, p_value, std_err = stats.linregress(stock_returns, mkt_returns)
beta, alpha
'''输出结果为 斜率、截距、相关系数、假设检验的p值、估算值的标准误差'''
通过回归估计出股票回报率对于市场变动的灵敏度$\beta$ ____$r_i = \alpha + \beta r_M$
然后通过：
$$E(R_i) = R_f + \beta_i [E(R_M) - R_f]$$
计算股票的期望回报率
同期实际回报率高于期望回报率，投资者在风险不变情况下获得更高的收益，则此证券价值被低估。
2.2套利定价模型
多元投资组合可以消除股票的非系统风险
套利定价理论 APT
假定证券的报酬率是基于多个系统风险因素的线性组合模型的到的（通货膨胀率、GDP增长率、实际利率或股息）
$$E(R_i) = \alpha_i + \beta_{i,1}F_1 + \beta_{i,2}F_2 + …+\beta_{i,j}F_j$$
要计算资产价格，必须计算所有的$\alpha_i$和$\beta$的值，因此需要在套利定价模型上进行多元线性回归。
import numpy as np
import statsmodels.api as sm
num_peiods = 9
all_values = np.array([np.random.random(8) for i in range(num_peiods)])
y_values = all_values[:, 0]
x_values = all_values[:, 1:]
x_values = sm.add_constant(x_values) #include the intercept
results = sm.OLS(y_values, x_values).fit()
results.summary()
2.4线性最优化
随着投资组合资产数量增加，限制条件也会增加，投资经理会在此限制下达到投资者的目标
线性最优化通过求解目标函数最小值或者最大值解决投资组合分配问题
开源的线性规划库——PuLP实现线性规划的单纯形法
https://github.com/coin-or/pulp
$$f(x, y) = 3x + 2y$$max
$$2x + y \leqslant 100$$
$$x + y \leqslant 80$$
$$x leqslant 40$$
$$x \geqslant 0, y \geqslant 0$$
import pulp
x = pulp.LpVariable('x', lowBound=0)
y = pulp.LpVariable('y', lowBound=0)
problem = pulp.LpProblem("A simple maximization objective", pulp.LpMaximize)
problem += 3*x + 2*y, 'The objective function'
problem += 2*x + y <= 100, '1st constraint'
problem += x + y <= 80, '2st constraint'
problem += x <= 40, '3rd constraint'
problem += x >= 0, '4th constraint'
problem += y >= 0, '5th constraint'
problem.solve()
for variable in problem.variables():
print(variable.name,'=', variable.varValue)
print(pulp.value(problem.objective))
整数规划
变量被限制为整数时成为整数规划——整数规划的特殊情况为二进制变量
整数规划模型经常用于实际问题建模
二进制条件下的整数规划模型
$$Minimize \sum_{i = x}^{z} IsOrder_i[variable \; cost_i \times quantity_i + fixed \; cost_i]$$
or
$$Minimize \sum_{i = x}^{z} variable \; cost_i \times quantity_i + fixed \; cost_i \times IsOrders$$
采用第一种实际计算了两个变量的乘法，但是这个库升级了，也可以直接计算。
建议采用第二种方法计算，它是真正的线性。
import pulp
dealers = ['X', 'Y', 'Z']
variable_costs = {'X': 500, 'Y': 450, 'Z': 450}
fixed_costs = {'X': 4000, 'Y': 2000, 'Z': 6000}
quantities = pulp.LpVariable.dicts('quantitity', dealers, lowBound=0, cat=pulp.LpInteger)
is_orders = pulp.LpVariable.dicts('orders', dealers, cat=pulp.LpBinary)
'''方法1'''
model = pulp.LpProblem('min', pulp.LpMinimize)
model += sum([(variable_costs[i] * quantities[i] + fixed_costs[i] * is_orders[i]) for i in dealers]), 'Minimize portfolio cost'
model += sum([quantities[i] for i in dealers]) == 150, 'total'
model += 30 <= quantities['X'] <= 100, 'boundary of total volume of X'
model += 30 <= quantities['Y'] <= 90, '..................Y'
model += 30 <= quantities['Z'] <= 70, '..................Z'
model.solve()
for variable in model.variables():
print(variable.name,'=', variable.varValue)
print(pulp.value(model.objective))
'''OUT'''
orders_X = 0.0
orders_Y = 0.0
orders_Z = 0.0
quantitity_X = 0.0
quantitity_Y = 90.0
quantitity_Z = 60.0
67500.0
##第二种计算方法
dealers = ['X', 'Y', 'Z']
variable_costs = {'X': 500, 'Y': 450, 'Z': 450}
fixed_costs = {'X': 4000, 'Y': 2000, 'Z': 6000}
quantities = pulp.LpVariable.dicts('quantitity', dealers, lowBound=0, cat=pulp.LpInteger)
is_orders = pulp.LpVariable.dicts('orders', dealers, cat=pulp.LpBinary)
model_2 = pulp.LpProblem('min', pulp.LpMinimize)
model_2 += sum([(variable_costs[i] * quantities[i] + fixed_costs[i] * is_orders[i]) for i in dealers]), 'Minimize portfolio cost'
model_2 += sum([quantities[i] for i in dealers]) == 150, 'total'
model_2 += is_orders['X'] * 30 <= quantities['X'] <= 100 * is_orders['X'], 'boundary of total volume of X'
model_2 += is_orders['Y'] * 30 <= quantities['Y'] <= 90 * is_orders['Y'], '..................Y'
model_2 += is_orders['Z'] * 30 <= quantities['Z'] <= 70 * is_orders['Z'], '..................Z'
model_2.solve()
for variable in model_2.variables():
print(variable.name,'=', variable.varValue)
print(pulp.value(model_2.objective))
'''OUT'''
orders_X = 0.0
orders_Y = 1.0
orders_Z = 1.0
quantitity_X = 0.0
quantitity_Y = 90.0
quantitity_Z = 60.0
75500.0
观察结果会发现，尽管第一种计算出了结果，但是它对应的x,y,z是否执行了并没有正确的给出。而且我通过excel同第二种结果一致
2.5 使用矩阵解线性方程组
$$2a+b+c=4$$
$$a+3b+2c=5$$
$$1a+0b+0c=6$$
$$\begin{bmatrix}
3 & 1 &1 \
1 & 3 &2 \
1 & 0 & 0
\end{bmatrix}
\times
\begin{bmatrix}
a\
b\
c
\end{bmatrix}
=
\begin{bmatrix}
4\
5\
6
\end{bmatrix}$$
$$Ax = B$$
$$x = A^{-1}B$$
import numpy as np
A = np.array([
[2, 1, 1],
[1, 3, 2],
[1, 0, 0]
])
B = np.array([4, 5, 6])
np.linalg.solve(A,B)
2.6 LU分解
$$A = LU$$
$L、U分别为下三角和上三角矩阵$
import scipy.linalg as linalg
A = np.array([
[2, 1, 1],
[1, 3, 2],
[1, 0, 0]
])
B = np.array([4, 5, 6])
LU = linalg.lu_factor(A)
x = linalg.lu_solve(LU, B)
x
P,L,U = scipy.linalg.lu(A)
print(P)
print(L)
print(U)
2.7 Cholesky分解
显著提高计算速度并降低对内存的要求
但是使用cholesky分解需要矩阵为埃尔米特矩阵（实值对称矩阵）且正定
$$A = LL^T$$
$L$为对角线为实数的下三角矩阵，$L^T$为$L$的共轭矩阵
import numpy as np
A = np.array([
[10, -1, 2, 0],
[-1, 11, -1, 3],
[2, -1, 10, -1],
[0, 3, -1, 8]
])
B = np.array([6, 25, -11, 15])
L = np.linalg.cholesky(A)
L
np.dot(L, L.T.conj())
$$L^T x = y = L^{-1} \times B$$
y = np.linalg.solve(L, B)
x = np.linalg.solve(L.T.conj(), y)
np.mat(A) * np.mat(x).T
2.8 QR分解
$$A = QR$$
$Q、R$分别是正交矩阵和上三角矩阵
$QR$算法是线性最小二乘问题的常见解法
正交矩阵乘以它的转置等于单位矩阵
正交矩阵的转置和它的逆向等
$$QRx = B$$
$$Rx = Q^{-1}B = Q^TB$$
A = np.array([
[2, 1, 1],
[1, 3, 2],
[1, 0, 0]
])
B = np.array([4, 5, 6])
Q, R = scipy.linalg.qr(A)
print(Q)
print(R)
y = np.dot(Q.T, B)
y
x = scipy.linalg.solve(R, y)
x
2.9 间接求解
Jacobi
$$A = D + R$$
$D$为对角矩阵
$$Ax = B$$
$$x_{n+1} = D^1(B - Bx_n)$$
def jacobi(A, B, n, tol=1e-20):
x = np.zeros_like(B)
for it_count in range(n):
x_new = np.zeros_like(x)
for i in range(A.shape[0]):
s1 = np.dot(A[i, :i], x[:i])
s2 = np.dot(A[i, i+1:], x[i+1:])
x_new[i] = (B[i] -s1 -s2) / A[i, i]
if np.allclose(x, x_new, tol):
break
x = x_new
return x
A = np.array([
[10.0, -1., 2., 0.],
[-1., 11., -1., 3.],
[2., -1., 10., -1.],
[0., 3., -1., 8.]
])
B = np.array([6., 25., -11., 15.])
n = 25
x = jacobi(A, B, n)
x
**程序中注意除法，输入的数据需要采用小数的形式。
这个问题出在python2与python3中的除法不一样中,包是按照2的习惯。
为了避免出现整数除法结果输出地板除结果，建议按照2的习惯，保持小数输入**
Gauss-Seidel
$$A = L + U$$
$L$为下三角矩阵,$U$为上三角矩阵
$$Ax = B$$
$$x_{n+1} = L^{-1}(B-Ux_n)$$
import numpy as np
def gauss(A, B, n, tol=1e-10):
L = np.tril(A)
U = A - L
L_inv = np.linalg.inv(L)
x = np.zeros_like(B)
for i in range(n):
Ux = np.dot(U, x)
x_new = np.dot(L_inv, B - Ux)
if np.allclose(x, x_new, tol):
break
x = x_new
return x
A = np.array([
[10.0, -1., 2., 0.],
[-1., 11., -1., 3.],
[2., -1., 10., -1.],
[0., 3., -1., 8.]
])
B = np.array([6., 25., -11., 15.])
n = 25
x = gauss(A, B, n)
x
博主个人能力有限，错误在所难免．
如发现错误请不要吝啬，发邮件给博主更正内容，在此提前鸣谢．
Email: JentChang@163.com （来信请注明文章标题，如果附带链接就更方便了）
你也可以在下方的留言板留下你宝贵的意见．建议邮件，它可以很快的被我发现，留言板的数据库我十年半月登录一次
