<jupyter_start><jupyter_text>Final Project<jupyter_code>import sys
import math
import numpy as np
import pandas as pd
import scipy.linalg as la
from statsmodels.tsa.ar_model import AR
from matplotlib import pyplot as plt
from sklearn.decomposition import fastica as sk_fastica
from sklearn.decomposition import PCA
from copy import deepcopy
import numba
from numba import jit<jupyter_output>/opt/conda/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
  from pandas.core import datetools<jupyter_text>Real Data Sets<jupyter_code>import pandas_datareader as pdr

data_source = 'morningstar'
#tickers = ['AAPL','AMD','CSCO','NVDA','ADBE','AMZN','EBAY','INTC','TSLA','FB','STX','TMUS','GOOGL','JD','MSFT','NFLX']
tickers = ['AAPL','AMD','CSCO','NVDA','ADBE','AMZN','EBAY','INTC']
start_date = '2016-01-01'
end_date = '2017-12-31'
data = pd.DataFrame()
for ticker in tickers:
    one_ticker = pdr.DataReader(ticker,data_source,start_date, end_date)
    one_ticker['Return'] = np.log(one_ticker['Close'].shift(1)/one_ticker['Close'])
    one_ticker = one_ticker.reset_index()[['Date','Return']].rename({'Return':ticker},axis=1).iloc[1:,:]
    if data.empty:
        data = one_ticker
    else:
        data = pd.merge(data, one_ticker, how='outer')
data.head()
data = data.fillna(0)
data.head()
pdata = pdr.DataReader(tickers,data_source,start_date, end_date)
pdata = pdata.reset_index()
plt.figure(figsize=(8,12))
i = 0
for ticker in tickers:
    plt.subplot(6, 2, i + 1)
    plt.plot(np.arange(len(pdata[pdata['Symbol']==ticker]['Close'])),pdata[pdata['Symbol']==ticker]['Close'])
    plt.title(ticker)
    i+=1
plt.tight_layout()
plt.savefig("prices.png")
plt.figure(figsize=(8,12))
for i in range(len(tickers)):
    plt.subplot(6, 2, i + 1)
    plt.plot(data.iloc[:,i+1])
    plt.title(data.columns.values[i+1])
plt.tight_layout()
plt.savefig("returns.png")
matrix_return = np.array(data.iloc[:,1:])
np.random.seed(123)<jupyter_output><empty_output><jupyter_text>FastICA Implementation<jupyter_code>#whiten
def whiten2(X):
    X = X - X.mean(-1)[:, None]
    A = np.dot(X, X.T)
    D, P = la.eigh(A)
    D = np.diag(D)
    S_inv = np.sqrt(np.linalg.pinv(D))
    V = np.dot(S_inv, P.T)
    return np.dot(V, X), V
def whiten(X, n_components):
    X = X - X.mean(-1)[:, None]
    u, d, _ = la.svd(X, full_matrices=False)
    V = (u / d).T[:n_components]
    X *= np.sqrt(X.shape[1])
    
    return np.dot(V, X), V
#three functions
def _logcosh(x):
    gx = np.tanh(x)
    g_x = 1- gx ** 2
    return gx, g_x.mean(-1)

def _exp(x):
    exp = np.exp(-(x ** 2) / 2)
    gx = x * exp
    g_x = (1 - x ** 2) * exp
    return gx, g_x.mean(axis=-1)

def _cube(x):
    return x ** 3, (3 * x ** 2).mean(axis=-1)

def decorrelation(W):
    s, u = la.eigh(np.dot(W, W.T))
    return np.dot(np.dot(u * (1. / np.sqrt(s)), u.T), W)

# fastICA
def fastICA(X, fun='logcosh',n_components=None, maxIter = 200, tol = 1e-04):
    X = X.T
    n, m = X.shape
    if n_components is None:
        n_components = n
    X1,V = whiten(X,n_components)
    p = float(m)
    if fun == 'logcosh':
        g = _logcosh
    elif fun == 'exp':
        g = _exp
    elif fun == 'cube':
        g = _cube
    
    #initialize w with normal distribution
    W = np.asarray(np.random.normal(size=(n_components,n_components)))

    # calculate w iteratively
  
    for ii in range(maxIter):
        gwtx, g_wtx = g(np.dot(W, X1))
        W1 = decorrelation(np.dot(gwtx, X1.T) / p - g_wtx[:, None] * W)
        lim = max(abs(abs(np.diag(np.dot(W1, W.T))) - 1))
        W = W1
        if lim < tol:
            break
    S = np.dot(np.dot(W, V), X).T
    return V,W,S<jupyter_output><empty_output><jupyter_text>Optimized Version<jupyter_code>@jit
def fastICA2(X, fun='logcosh',n_components=None, maxIter = 200, tol = 1e-04):
    X = X.T
    n, m = X.shape
    if n_components is None:
        n_components = n
    X1,V = whiten(X,n_components)
    p = float(m)
    if fun == 'logcosh':
        g = _logcosh
    elif fun == 'exp':
        g = _exp
    elif fun == 'cube':
        g = _cube
    
    #initialize w with normal distribution
    W = np.asarray(np.random.normal(size=(n_components,n_components)))

    # calculate w iteratively
  
    for ii in range(maxIter):
        gwtx, g_wtx = g(np.dot(W, X1))
        W1 = decorrelation(np.dot(gwtx, X1.T) / p - g_wtx[:, None] * W)
        lim = max(abs(abs(np.diag(np.dot(W1, W.T))) - 1))
        W = W1
        if lim < tol:
            break
    S = np.dot(np.dot(W, V), X).T
    return V,W,S<jupyter_output><empty_output><jupyter_text>Application to simulated data Generate simulated data<jupyter_code>class simulated_data():
    def __init__(self, signalRange, timeRange):
        self.signalRange = signalRange
        self.timeRange = timeRange
        self.x = np.arange(0, self.timeRange, 0.1)
        self.pointNumber = self.x.shape[0]

    def produceSin(self, period=100):
        y = self.signalRange * np.sin(self.x / period * 2 * np.pi)
        return y

    def produceTri(self, period=20):
        lastPoint = period * 10 - 1
        if lastPoint >= self.pointNumber:
            raise ValueError('At least one period!')
        delta = ((-1 - 1) * self.signalRange) / float(self.x[lastPoint] - self.x[0])
        y = (self.x[:lastPoint+1] - self.x[0]) * delta + self.signalRange
        y = np.tile(y, self.pointNumber // y.shape[0])[:self.pointNumber]
        return -y

    def mixSignal(self,majorSignal, *noises):
        mixSig = deepcopy(majorSignal)
        noiseRange = 0.5
        for noise in noises:
            mixSig += noiseRange * noise
        return mixSig
sdata = simulated_data(2, 100)
Sin = sdata.produceSin(15)
Tri = sdata.produceTri(25)
totalSig = [Sin, Tri]
mixSig = []
for i, majorSig in enumerate(totalSig):
    curSig = sdata.mixSignal(majorSig, *(totalSig[:i] + totalSig[i+1:]))
    mixSig.append(curSig)
mixSig = np.asarray(mixSig)<jupyter_output><empty_output><jupyter_text>The Source Signals<jupyter_code>signals = np.c_[Sin,Tri]
plt.figure(figsize=(10,3))
for i in range(2):
    plt.subplot(2, 1, i + 1)
    plt.plot(np.arange(0, 100, 0.1),signals[:,i])
plt.tight_layout() 
plt.savefig('source.png')<jupyter_output><empty_output><jupyter_text>The observed mixtures of the source signal<jupyter_code>plt.figure(figsize=(10,3))
for i in range(2):
    plt.subplot(2, 1, i + 1)
    plt.plot(np.arange(0, 100, 0.1),mixSig[i])
plt.tight_layout()
plt.savefig('mixed.png')<jupyter_output><empty_output><jupyter_text>Apply ICA to the mixed signals<jupyter_code>#do ica by using our implementation
V,W,S = fastICA(mixSig.T, fun='logcosh',maxIter = 800)
m,n = S.T.shape
#print(S)
plt.figure(figsize=(10,3))
for i in range(2):
    plt.subplot(2, 1, i + 1)
    plt.plot(np.arange(0,100,0.1), S.T[i,:])
plt.tight_layout()
plt.savefig("ICs.png")<jupyter_output><empty_output><jupyter_text>Using sklearn package<jupyter_code>V,W,S = sk_fastica(mixSig.T, fun='logcosh')
m,n = S.T.shape
#print(S)
plt.figure(figsize=(10,3))
for i in range(2):
    plt.subplot(2, 1, i + 1)
    plt.plot(np.arange(0,100,0.1), S.T[i,:])
plt.tight_layout()
plt.savefig("ICs_sklearn.png")<jupyter_output><empty_output><jupyter_text>Application to real data Do ICA without assigning value to n_components<jupyter_code>V,W,S = sk_fastica(matrix_return)
plt.figure(figsize=(12,12))
m,n = S.T.shape
for i in range(m):
    plt.subplot(m, 2, i + 1)
    plt.plot(np.arange(n), S.T[i])
plt.tight_layout()
plt.savefig("skpackage_all.png")
#do ica by using our own package
V,W,S = fastICA(matrix_return, fun='logcosh')
#S = np.dot(np.dot(W, V), matrix_return.T)
m,n = S.T.shape
#print(S)
plt.figure(figsize=(12,12))
for i in range(m):
    plt.subplot(m, 2, i + 1)
    plt.plot(np.arange(n), S.T[i])
plt.tight_layout()
plt.savefig("ownpackage_all.png")
x_hat = np.linalg.pinv(np.dot(W,V)) @ S.T
plt.plot(x_hat[1]-0.4,color = 'black')
plt.plot(matrix_return.T[1],color='blue')
plt.yticks([])
plt.legend(['Recovered','Original'])
#np.mean((x_hat[1]*100-matrix_return.T[1,:]*100)**2)
plt.title('AMD')
plt.tight_layout()
plt.savefig("compare.png")<jupyter_output><empty_output><jupyter_text>Let n_components = 4<jupyter_code>V,W,S = sk_fastica(matrix_return,n_components=4)
plt.figure(figsize=(12,12))
m,n = S.T.shape
for i in range(m):
    plt.subplot(m, 2, i + 1)
    plt.plot(np.arange(n), S.T[i])
plt.tight_layout()
plt.savefig("skpackage.png")
#do ica by using our own package
#mw, V = whiten(matrix_return.T,n_components=4)
V,W,S = fastICA(matrix_return, fun='logcosh',n_components=4)
#S = np.dot(np.dot(W, V), matrix_return.T)
m,n = S.T.shape
#print(S)
plt.figure(figsize=(12,12))
for i in range(m):
    plt.subplot(m, 2, i + 1)
    plt.plot(np.arange(n), S.T[i])
plt.tight_layout()
plt.savefig("ownpackage.png")<jupyter_output><empty_output><jupyter_text>Time series prediction with Independent Component Analysis<jupyter_code>def predict_ICA(S,q):
    n, m = S.shape
    s_pred = np.zeros((n,m-q))
    for i in range(n):
        model = AR(S[i,:])
        model_fit = model.fit(q)
        # make predictions
        predictions = model_fit.predict()
        s_pred[i,:] = predictions
    return s_pred

def text_mse(x_hat, x_real, n, q):
    mse = np.zeros((2,n))
    for i in range(n):
        #mse_ica
        mse[0,i] = np.sum((x_hat[i] - x_real[i,q:])**2)/(len(x_hat[i]))
        #mse_ar
        armodel = AR(x_real[i,:]).fit(q)
        mse[1,i] = np.sum((armodel.predict() - x_real[i,q:])**2)/(len(armodel.predict()))
    return mse

def compare_ica_ar(S, x_real, A, qs):
    n = x_real.shape[0]
    df = pd.DataFrame()
    for q in qs:
        sp = predict_ICA(S,q)
        x_hat = np.dot(A, sp)
        mse = text_mse(x_hat, x_real, n, q)*10000
        df = df.append(pd.DataFrame(mse))
        #df.columns = tickers
        #print("AR",q)
        #print("--"*20)
    df.columns = tickers
    df.to_csv('table.csv')
A = np.linalg.pinv(np.dot(W,V))
qs = [1,2,3,4,5,20]
compare_ica_ar(S.T, matrix_return.T, A, qs)<jupyter_output><empty_output><jupyter_text>PCA vs ICA<jupyter_code>pca = PCA(n_components=4)
pca_data = pca.fit_transform(matrix_return)
m = pca_data.shape[1]
plt.figure(figsize=(12,12))
for i in range(m):
    plt.subplot(m, 2, i + 1)
    plt.plot(pca_data[:,i],color='black')
    if i==0:
        j=0
    if i==1:
        j=3
    if i==2:
        j=1
    if i==3:
        j=2
    plt.plot(S.T[j]-0.4,color='blue')
    plt.legend(['PCA','ICA'])
plt.tight_layout()
plt.savefig("pca.png")<jupyter_output><empty_output><jupyter_text>Other Algorithms<jupyter_code>%%timeit
V,W,S = sk_fastica(matrix_return,n_components=4, algorithm = 'parallel')
%%timeit
V,W,S = sk_fastica(matrix_return,n_components=4, algorithm='deflation')
%%timeit
V,W,S = fastICA(matrix_return,n_components=4)
%%timeit
V,W,S = fastICA2(matrix_return,n_components=4)<jupyter_output>3.66 ms ± 235 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
{'fraction_non_alphanumeric': 0.1706785357434826, 'fraction_numerical': 0.02425288400399673, 'mean_word_length': 3.67714528462192, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '357996', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4738, 'n_tokens_neox': 4517, 'n_words': 864, 'path': 'Final Project/Final Project-Zhechang Yang and Xi Chen.ipynb', 'repo_name': 'ZhechangYang/STA663', 'size': '11009'}