# 读入"指数行情序列daily.xlsx"
stockinfo=pd.read_excel("指数行情序列daily.xlsx",index_col = 0, skiprows = 1,na_values=['0.0000'])
stockinfo.head()

# '沪深300','上证国债','人民币实际指数','Wind商品指数' 的月对数收益率
range_month = pd.date_range('1990-12-19', '2017-10-18', freq='BM')
stockinfo_month = DataFrame(stockinfo,index = range_month, columns=['沪深300','上证国债','人民币实际指数','Wind商品指数'])
stockinfo_month = stockinfo_month.fillna(method='ffill')

logp_month=stockinfo_month.applymap(log)
LogReturnMonth=logp_month-logp_month.shift(1)

LogReturnMonth.head()

#投资组合分析
LogReturnMonth.describe()

Mu=LogReturnMonth.mean()
Sigma=LogReturnMonth.cov()
Sigma

Sigma=np.matrix(Sigma)
N=max(Mu.shape)
Mu=np.matrix(Mu).T
onesv=np.matrix(np.ones((N,1)))
SI=Sigma.I
Mu.T*SI*onesv
A=Mu.T*SI*onesv; A=A[0,0]
B=Mu.T*SI*Mu; B=B[0,0]
C=onesv.T*SI*onesv; C=C[0,0]
D=B*C-A**2
g=SI*onesv*B/D-SI*Mu*A/D
h=SI*Mu*C/D-SI*onesv*A/D
mmin=np.min(Mu)
mmax=np.max(Mu)
meanv=np.linspace(mmin,mmax,100)
wv=list(map(lambda x: g+x*h,meanv))
sigmav=list(map(lambda x: np.sqrt(x.T*Sigma*x)[0,0],wv))
mumin=-g.T*Sigma*h/(h.T*Sigma*h)
sigmamin=np.sqrt(g.T*Sigma*g-(g.T*Sigma*h)**2/(h.T*Sigma*h))
muf=0.001    #设无风险收益率为 0.001
wt=SI*(Mu-muf)
wt=wt/(onesv.T*wt)
mut=wt.T*Mu
sigmat=np.sqrt(wt.T*Sigma*wt)

%matplotlib inline
Sig=LogReturnMonth.std()
plt.plot(sigmav,meanv)
plt.plot(Sig,Mu,"bo")
plt.plot(sigmamin,mumin,"ro")


plt.plot([0,sigmat.item(0)],[muf,mut.item(0)],"yo--")


plt.xlim([-0.005,0.09]) 
plt.ylim([-0.001,0.013])

LogReturnMonth.describe()

LogReturnMonth1 = DataFrame(LogReturnMonth,columns = ['沪深300','人民币实际指数','Wind商品指数'])
Mu=LogReturnMonth1.mean()
Sigma=LogReturnMonth1.cov()
Sigma

Sigma=np.matrix(Sigma)
N=max(Mu.shape)
Mu=np.matrix(Mu).T
onesv=np.matrix(np.ones((N,1)))
SI=Sigma.I

Mu.T*SI*onesv

A=Mu.T*SI*onesv; A=A[0,0]
B=Mu.T*SI*Mu; B=B[0,0]
C=onesv.T*SI*onesv; C=C[0,0]
D=B*C-A**2
g=SI*onesv*B/D-SI*Mu*A/D
h=SI*Mu*C/D-SI*onesv*A/D

mmin=np.min(Mu)
mmax=np.max(Mu)
meanv=np.linspace(mmin,mmax,100)
wv=list(map(lambda x: g+x*h,meanv))
sigmav=list(map(lambda x: np.sqrt(x.T*Sigma*x)[0,0],wv))

mumin=-g.T*Sigma*h/(h.T*Sigma*h)
sigmamin=np.sqrt(g.T*Sigma*g-(g.T*Sigma*h)**2/(h.T*Sigma*h))

muf=0.002668           #无风险收益率
wt=SI*(Mu-muf)
wt=wt/(onesv.T*wt)
mut=wt.T*Mu
sigmat=np.sqrt(wt.T*Sigma*wt)

Sig=LogReturnMonth1.std()
plt.plot(sigmav,meanv)
plt.plot(Sig,Mu,"bo")
plt.plot(sigmamin,mumin,"ro")


plt.plot([0,sigmat.item(0)],[muf,mut.item(0)],"yo--")


plt.xlim([-0.005,0.10]) 
plt.ylim([0,0.013])

weightm=np.array([[x.item(0),x.item(1),x.item(2)] for x in wv])
for i in range(3):
    plt.plot(meanv,weightm[:,i],label=i+1)
plt.legend(loc="best")
plt.title("weighs for 3 assets as function of $\mu_p$")

#先考虑第一组的数据
LogReturnMonth_1.describe()

##########  第一部分没有 ‘人民币实际指数 ’的数据。所以只考虑 '沪深300'和'Wind商品指数'的组合
LogReturnMonth_1 = DataFrame(LogReturnMonth_1,columns = ['沪深300','Wind商品指数'])
Sigma1=LogReturnMonth_1.cov()
Sigma1

import sympy as sy
w=sy.symbols('w',positive = True, real = True)
ExpRet=mu2+(mu1-mu2)*w
Varf=sy.simplify(sig1**2*w**2+sig2**2*(1-w)**2+2*w*(1-w)*rho12*sig1*sig2)
StdRet=sy.sqrt(Varf)
#init_printing()
ExpRet,sy.simplify(StdRet)

meanf= sy.lambdify(w, ExpRet, "numpy")
varf=sy.lambdify(w,Varf,"numpy")

#%matplotlib qt
%matplotlib inline

fig=plt.figure(figsize=(8,6)) 

plt.text(sig1+0.002,mu1-0.001,"$R_1$")
plt.text(sig2+0.002,mu2-0.001,"$R_2$")

wmin=float(-Varf.coeff(w,1)/(2*Varf.coeff(w,2)))
mx=np.sqrt(varf(wmin));my=meanf(wmin)
plt.plot([mx,0],[my,muf],"bo")
plt.text(mx+0.005,my,"M (minimum variance portfolio)",fontsize=15)
plt.text(0,muf+0.005,"F",fontsize=15)

wv=np.linspace(0,wmin,100)
stdv=np.sqrt(varf(wv))
meanv=meanf(wv)
plt.plot(stdv,meanv,"r-")  
plt.annotate("$w\in$"+str([0,wmin]), size=15,xy=(stdv[50]+0.002,meanv[50]),
             xytext=(stdv[50]+0.05,meanv[50]),
            arrowprops=dict(facecolor="red"))


wv=np.linspace(wmin,1,100)
stdv=np.sqrt(varf(wv))
meanv=meanf(wv)
plt.plot(stdv,meanv,"b-") 
plt.annotate("efficient frontier\n"+"$w\in$"+str([wmin,1]), size=15,xy=(stdv[80]+0.002,meanv[80]),
             xytext=(stdv[80]+0.05,meanv[80]),
            arrowprops=dict(facecolor="blue"))

plt.xlim([-0.02,max(sig1,sig2)*1.5])
plt.ylim([-0.008,0.002])
plt.plot([0,mx],[muf,my],"b--")
plt.plot([0,sig2],[muf,mu2],"r-.")

# tangency portfolio
v1=mu1-muf;v2=mu2-muf;
Wt=(v1*sig2**2-v2*rho12*sig1*sig2)/(v1*sig2**2+v2*sig1**2-(v1+v2)*rho12*sig1*sig2)
Tx=np.sqrt(varf(Wt));Ty=meanf(Wt)
plt.plot([Tx],[Ty],"bo")
plt.plot([Tx,0],[Ty,muf],"b--")
plt.text(Tx+0.005,Ty,"T (tangency portfolio)",fontsize=15)

#然后虑第二组的数据
LogReturnMonth_2.describe()

#无风险资产收益率为 0.003452
LogReturnMonth_2 = DataFrame(LogReturnMonth_2,columns = ['沪深300','人民币实际指数','Wind商品指数'])
Mu=LogReturnMonth_2.mean()
Sigma=LogReturnMonth_2.cov()
Sigma

Sigma=np.matrix(Sigma)
N=max(Mu.shape)
Mu=np.matrix(Mu).T
onesv=np.matrix(np.ones((N,1)))
SI=Sigma.I

Mu.T*SI*onesv

A=Mu.T*SI*onesv; A=A[0,0]
B=Mu.T*SI*Mu; B=B[0,0]
C=onesv.T*SI*onesv; C=C[0,0]
D=B*C-A**2
g=SI*onesv*B/D-SI*Mu*A/D
h=SI*Mu*C/D-SI*onesv*A/D

mmin=np.min(Mu)
mmax=np.max(Mu)
meanv=np.linspace(mmin,mmax,100)
wv=list(map(lambda x: g+x*h,meanv))
sigmav=list(map(lambda x: np.sqrt(x.T*Sigma*x)[0,0],wv))

mumin=-g.T*Sigma*h/(h.T*Sigma*h)
sigmamin=np.sqrt(g.T*Sigma*g-(g.T*Sigma*h)**2/(h.T*Sigma*h))

muf=0.003452          #无风险收益率
wt=SI*(Mu-muf)
wt=wt/(onesv.T*wt)
mut=wt.T*Mu
sigmat=np.sqrt(wt.T*Sigma*wt)

Sig=LogReturnMonth_2.std()
plt.plot(sigmav,meanv)
plt.plot(Sig,Mu,"bo")
plt.plot(sigmamin,mumin,"ro")


plt.plot([0,sigmat.item(0)],[muf,mut.item(0)],"yo--")


plt.xlim([-0.005,0.10]) 
plt.ylim([0,0.013])

weightm=np.array([[x.item(0),x.item(1),x.item(2)] for x in wv])
for i in range(3):
    plt.plot(meanv,weightm[:,i],label=i+1)
plt.legend(loc="best")
plt.title("weighs for 3 assets as function of $\mu_p$")

stocks=pd.read_csv("stockclose.csv",index_col = 0,parse_dates = True)
stocks

range_week = pd.date_range('2014-01-01', '2018-12-31', freq='W-FRI')
stocks = DataFrame(stocks, index = range_week, columns = ['000001.SZ','000002.SZ','000004.SZ','000005.SZ','000006.SZ',
                                                          '000007.SZ','000008.SZ','000009.SZ','000010.SZ','000011.SZ'])
                   
logp_week=stocks.applymap(log)
LogReturnWeek=logp_week-logp_week.shift(1)
LogReturnWeek = LogReturnWeek.fillna(method = 'ffill')
LogReturnWeek

# 随机取两支股票
from random import randint
columns = ['000001.SZ','000002.SZ','000004.SZ','000005.SZ','000006.SZ','000007.SZ','000008.SZ','000009.SZ','000010.SZ','000011.SZ']
x = randint(0,len(columns)-1)
random_stock1 = columns[x]
columns.pop(x)
y = randint(0,len(columns)-1)
random_stock2 = columns[y]
columns.pop(y)

New_Logreturn = DataFrame(LogReturnWeek,columns = [random_stock1,random_stock2])
New_Logreturn.head()

New_Logreturn_describe = New_Logreturn.describe()
New_Logreturn_describe

New_Logreturn_cov = New_Logreturn.cov()
New_Logreturn_cov

mu1=New_Logreturn_describe.iloc[1][0]
mu2=New_Logreturn_describe.iloc[1][1]
sig1=New_Logreturn_describe.iloc[2][0]
sig2=New_Logreturn_describe.iloc[2][1]
rho12=New_Logreturn_cov.iloc[0][1]
ExpRet=mu2+(mu1-mu2)*w
Varf=sy.simplify(sig1**2*w**2+sig2**2*(1-w)**2+2*w*(1-w)*rho12*sig1*sig2)
StdRet=sy.sqrt(Varf)
ExpRet,sy.simplify(StdRet)

meanf= sy.lambdify(w, ExpRet, "numpy")
varf=sy.lambdify(w,Varf,"numpy")

%matplotlib inline

fig=plt.figure(figsize=(8,6)) 

plt.text(sig1+0.002,mu1-0.001,"$R_1$")
plt.text(sig2+0.002,mu2-0.001,"$R_2$")

wmin=float(-Varf.coeff(w,1)/(2*Varf.coeff(w,2)))
mx=sy.sqrt(varf(wmin));my=meanf(wmin)
plt.plot(mx,my,"bo")
plt.text(mx+0.005,my,"minimum variance portfolio")


wv=np.linspace(0,wmin,100)
stdv=np.sqrt(varf(wv))
meanv=meanf(wv)
plt.plot(stdv,meanv,"r-")  
plt.annotate("$w\in$"+str([0,wmin]), size=20,xy=(stdv[50]+0.002,meanv[50]),
             xytext=(stdv[50]+0.05,meanv[50]),
            arrowprops=dict(facecolor="red"))

plt.xlim([0,0.09])
plt.ylim([-0.006,0.01])
plt.xlabel(r"$\sigma_r$")
plt.ylabel(r"$E(R)$")
wv=np.linspace(wmin,1,100)
stdv=np.sqrt(varf(wv))
meanv=meanf(wv)
plt.plot(stdv,meanv,"b-") 
plt.annotate("efficient frontier\n"+"$w\in$"+str([wmin,1]), size=20,xy=(stdv[50]+0.002,meanv[50]),
             xytext=(stdv[50]+0.05,meanv[50]),
            arrowprops=dict(facecolor="blue"))

print("由于随机取2只股票，不知道图出现在哪儿，所以把坐标系设大一点")

# 为了方便按columns里的股票循序依次添加
columns