815主成分分析笔记

开始之前，让我们一起说：“谢谢Katerina”

好的，本期笔记是tutorial 2的，代码内容是question 1，question 2自己看答案吧，毕竟它也没有代码可以写。
首先是第一部分，是一个矩阵分解的函数，说是函数呢，矩阵分解是用numpy做的，这个函数的输出也不清楚在做什么，甚至后面也没有用到这个函数，无语

# -*- coding: utf-8 -*-
"""
Created on Fri Feb  4 13:19:33 2022

@author: Pamplemousse
"""
# In[1]:

import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from matplotlib import pyplot as plt
sns.set()

def factors_XX_SVD(x,k, stand=1):#stand=1，必须标准化的意思
    if stand:
        sc=StandardScaler(copy=True, with_mean=True, with_std=False)
        x=sc.fit_transform(x)#标准化
    
    U,S,V=np.linalg.svg(x,full_matrices=False) #矩阵分解
    #如果x是T*N，那么U是T*T,S是T*N,V是N*N，x=U*S*V
    T,N=x.shape#如果X是a*b，T=a,N=b
    S=np.diag(S)  #对角线元素
    
    loadings=(1/np.sqrt(N))*U[:,:k]@S[:k,:k]#U的前k列和S的前k个对角元素矩阵相乘，除以sqrt(N)
    #注意，前面有说N是x的列数
    
    factors=np.sqrt(N)*V[:k,:]#V的前k行除以sqrt(N)
    
    recon=U[:,:k]@S[:k,:k]@V[:k,:]#前k个主成分的拟合结果
    
    return factors, loadings, recon
#就是说这个函数不知道在这里干嘛的，后面也没有用，这些输出是拿来干嘛的咱也不知道
#离谱 （对，我就是为了说一句离谱才把它复制上的）

接下来是question1的第一个小点，主要是对股票收益率做主成分分析。虽然老师把它分成了如下这么多块，但是下面的代码（连着上面的也可以）一起跑了也没什么问题。

# In[2]:
    
df=pd.read_csv('intradayreturns.csv')
#这个.csv文件是课程网页上下载的，是使用pandas读取的，不一定要命名成df，但这个数据的类型是dataframe
n,t=df.shape#获取df的维度


# In[3]:

hld=df.to_numpy()#将dataframe转成numpy的数据类型
x=np.zeros([n,t])#设置一个大小为n*t的0矩阵，用于之后存放股票收益率数据

for i in range(0,n):
    x[i,:]=(hld[i,]-np.mean(hld[i,]))/np.std(hld[i,])  #手动标准化



# In[4]:
    
U,S,Vt=np.linalg.svd(x,full_matrices=False,compute_uv=True)#使用numpy的函数进行矩阵分解，具体分解原理自行百度。
s=np.diag(S)


# In[5]:

Rmax=np.min([n,t])   #满秩时秩的大小
plt.plot(np.asarray(list(range(Rmax))),S[:Rmax],'*');
plt.title('Scree Plot')
plt.xlabel('PC number')
plt.ylabel('Singular Value')

ratio=s[0,0]/s[1,1]
print('Ratio of the first over the second principal components is ', ratio)

得到主成分图像：

其中， Ratio of the first over the second principal components is 2.2899386906267933

下面是第二个小点，画出前x个主成分对数据总的解释效力。

# In[6]:

# Question B
fa=np.zeros([Rmax+1,1])#fa是一个大小为（秩+1）*1的0矩阵
for pick in range(1,Rmax+1):#i遍历从1到Rmax
    x2=U[:,0:pick]@s[:pick,:pick]@Vt[:pick,:] #x2是拟合的数据
    res=x-x2 #实际-拟合=残差
    fa[pick]=1-(np.power(res.flatten(),2).sum())/(np.power(x,2).sum())
    #1-残差平方和/总平方和，即解释平方和
    #因为x经过了标准化，是0均值的，所以直接平方就是总平方和

fa=fa[1:] #fa[0]是0
plt.plot(np.asarray(list(range(Rmax))),fa,'*')
plt.title('Fraction Explained')
plt.xlabel('PC number')
plt.ylabel('Fraction Explained')

fa[0]#输出第一个主成分的解释性大小

# In[7]:
#或者更规范一点，这样输出第一个主成分的解释平方和：
fraction_explained_all_R1=fa[0]
print(fraction_explained_all_R1)

得到前N个主成分的解释力图像：

其中第一个主成分的解释平方和为fa[0]，输出结果是0.12083467

最后是第三个小问，绘制第1、3、5个主成分对每一支股票的解释性。

# In[8]:

#Question C
R=[1,3,5]
q=0
fa=np.zeros([3,13424])
#x.shape是（13424,245）

for pick in R:
    x2=U[:,pick-1:pick]@s[pick-1:pick,pick-1:pick]@Vt[pick-1:pick,:]
    #老师给的代码这里是前pick个主成分，我觉得应该改成这样
    #但是，老师的结果跟他的解释是完全对不上的，莫名其妙
    res=x-x2
    #这里求残差
    
    for ii in range(0,n):
        fa[q,ii]=1-np.power(res[ii,:],2).sum()/np.power(x[ii,:],2).sum()
        #对每一个pick的第i列求解释平方和，一共求3*13424个解释平方和
    
    fa[q,:]=np.sort(fa[q,:])#对3个解释平方和序列分别从小到大排序
    q+=1


# In[9]:
    
plt.plot(fa[0,:],'b-') #第1个主成分-蓝色
plt.plot(fa[1,:],'g-') #第3个主成分-绿色
plt.plot(fa[2,:],'r-') #第5个主成分-红色
plt.title('fraction explained for individual stock, R=1,3,5')
plt.xlabel('stock number (sorted by fraction explained)')
plt.ylabel('fraction explained for that stock')

可以想象，总的来说，蓝色的线在最上面，其次是绿色红色，因为排在前面的主成分解释性更高嘛。得到如下图像：（老师真的应该给图片搞个图例，因为他没有搞，所以Katerina乱讲，很烦。但是我也懒得加个图例 -_-!）

可以看到绿色的线确实xue微比红色的高哈，因为主成分是各个变量合成的，所以靠后的主成分对单独的变量最高的解释平方和要低于前面的主成分。（反正这个才更符合第pick个主成分对单独变量的解释平方和图像，我个人是这样认为的。）

如果发现什么缺漏或者bug请联系我。
现在，下班！