<jupyter_start><jupyter_text>KALMAN FILTER<jupyter_code>import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

def Kalman_Filter(Y):
     S = Y.shape[0]
     S = S + 1
     Z = param0[0]
     T = param0[1]
     H = param0[2]
     Q = param0[3]
     u_predict = np.zeros(S)
     u_update = np.zeros(S)
     P_predict = np.zeros(S)
     P_update = np.zeros(S)
     v = np.zeros(S)
     F = np.zeros(S)
     KF_Dens = np.zeros(S)
     for s in range(1,S):
      if s == 1: 
        P_update[s] = 1000
        P_predict[s] =  T*P_update[1]*np.transpose(T)+Q    
      else: 
        F[s] = Z*P_predict[s-1]*np.transpose(Z)+H 
        v[s]=Y[s-1]-Z*u_predict[s-1]   
        u_update[s] = u_predict[s-1]+P_predict[s-1]*np.transpose(Z)*(1/F[s])*v[s]
        u_predict[s] = T*u_update[s]; 
        P_update[s] = P_predict[s-1]-P_predict[s-1]*np.transpose(Z)*(1/F[s])*Z*P_predict[s-1];
        P_predict[s] = T*P_update[s]*np.transpose(T)+Q
        KF_Dens[s] = (1/2)*np.log(2*np.pi)+(1/2)*np.log(abs(F[s]))+(1/2)*np.transpose(v[s])*(1/F[s])*v[s]     
        Likelihood = np.sum(KF_Dens[1:-1])  
        return Likelihood


def Kalman_Smoother(params, Y):
 S = Y.shape[0]
 S = S + 1
 "Initialize Params:"
 Z = params[0]
 T = params[1]
 H = params[2]
 Q = params[3]
 
 "Kalman Filter Starts:"
 u_predict = np.zeros(S)
 u_update = np.zeros(S)
 P_predict = np.zeros(S)
 P_update = np.zeros(S)
 v = np.zeros(S)
 F = np.zeros(S)
 for s in range(1,S):
   if s == 1: 
    P_update[s] = 1000
    P_predict[s] =  T*P_update[1]*np.transpose(T)+Q    
   else: 
    F[s] = Z*P_predict[s-1]*np.transpose(Z)+H 
    v[s]=Y[s-1]-Z*u_predict[s-1]   
    u_update[s] = u_predict[s-1]+P_predict[s-1]*np.transpose(Z)*(1/F[s])*v[s]
    u_predict[s] = T*u_update[s]; 
    P_update[s] = P_predict[s-1]-P_predict[s-1]*np.transpose(Z)*(1/F[s])*Z*P_predict[s-1];
    P_predict[s] = T*P_update[s]*np.transpose(T)+Q
    
    u_smooth = np.zeros(S)
    P_smooth = np.zeros(S)
    u_smooth[S-1] = u_update[S-1]
    P_smooth[S-1] = P_update[S-1]    
 for  t in range(S-1,0,-1):
        u_smooth[t-1] = u_update[t] + P_update[t]*np.transpose(T)/P_predict[t]*(u_smooth[t]-T*u_update[t])
        P_smooth[t-1] = P_update[t] + P_update[t]*np.transpose(T)/P_predict[t]*(P_smooth[t]-P_predict[t])/P_predict[t]*T*P_update[t]
 u_smooth = u_smooth[0:-1]
 return u_smooth
T = 861
Y = openmean
param0 = np.array([0.5, 0.5, openmean.var()*100, 100*openmean.var()])
param_star = minimize(Kalman_Filter, param0, method='BFGS', options={'xtol': 1e-8, 'disp': True})
Y_update = Kalman_Smoother(param_star.x, Y)
timevec = np.linspace(1,T,T)
plt.plot(timevec[3:-1], Y_update[3:-1],'r',timevec, Y,'b:')
def rmseCalc(Y,Y_hat):
    rmse = np.sqrt(np.mean((Y_update-Y)**2))
    print(rmse)
rmseCalc(Y,Y_update)<jupyter_output>4.400720734460765<jupyter_text>GARCH<jupyter_code>T = 861;
mu = openmean.mean();
sig  = openmean.var();
Y =openmean
def GARCH(Y):
    mu = param0[0]
    omega = np.exp(param0[1])
    alpha = (1-np.exp(-param0[2])) / (1+np.exp(-param0[2]))
    beta = param0[3]
    T = Y.shape[0]
    GARCH_Dens = np.zeros(T)
    sigma2 = np.zeros(T)
    F = np.zeros(T)
    v = np.zeros(T)
    
    for t in range(1,T):
        sigma2[t] = omega+alpha*((Y[t-1]-mu)**2)+beta*(sigma2[t-1]);
        F[t] = Y[t] - mu-np.sqrt(sigma2[t])*np.random.normal(0,1,1);
        v[t] = sigma2[t];
        GARCH_Dens[t] = (1/2)*np.log(2*np.pi)+(1/2)*np.log(v[t])+(1/2)*(F[t]/v[t])
    
    Likelihood = np.sum(GARCH_Dens[1:-1])
    return Likelihood
def GARCH_PROD(params, Y0, T):
    mu = params[0]
    omega = np.exp(params[1])
    alpha = (1-np.exp(-params[2])) / (1+np.exp(-params[2]))
    beta = params[3]
    Y = np.zeros(T)
    sigma2 = np.zeros(T)
    Y[0] = Y0
    sigma2[0] = omega/(1-alpha)
    for t in range (1,T):
        sigma2[t] = omega+alpha*((Y[t-1]-mu)**2)+beta*(sigma2[t-1]);
        Y[t] = mu-np.sqrt(sigma2[t])*np.random.normal(0,1,1);
    return Y  
param_star = minimize(GARCH, param0, method = 'BFGS',options = {'xtol':1e-8,'disp':True})
#param_star = minimize(Kalman_Filter, param0, method='BFGS', options={'xtol':1e-8, 'disp': True})
param0 = np.array([mu,2.5,0.3,0.5])
param_star = minimize(GARCH, param0, method = 'BFGS',options = {'xtol':1e-8,'disp':True})
Y_GARCH = GARCH_PROD(param_star.x,Y[0],T)
#Y_GARCH = GARCH_PROD(param_star.x,Y[0],T)
timevec = np.linspace(1,T,T)
plt.plot(timevec, Y ,'b', timevec[5:-1], Y_GARCH[5:-1],'r:')
rmse = np.sqrt(np.mean((Y_GARCH-Y)**2))
print(rmse)<jupyter_output>26.812712148714063<jupyter_text>2 weeks data fitting for Garch and Garch-T<jupyter_code>fb=facebook[350:379]
fb.head()
fb.shape

openmean1= (fb['Open']+ fb['Close'])/2
#openmean

def GARCH(Y):
    mu = param0[0]
    omega = np.exp(param0[1])
    alpha = (1-np.exp(-param0[2])) / (1+np.exp(-param0[2]))
    beta = param0[3]
    T = Y.shape[0]
    GARCH_Dens = np.zeros(T)
    sigma2 = np.zeros(T)
    F = np.zeros(T)
    v = np.zeros(T)
    
    for t in range(1,T):
        sigma2[t] = omega+alpha*((Y[t-1]-mu)**2)+beta*(sigma2[t-1]);
        F[t] = Y[t] - mu-np.sqrt(sigma2[t])*np.random.normal(0,1,1);
        v[t] = sigma2[t];
        GARCH_Dens[t] = (1/2)*np.log(2*np.pi)+(1/2)*np.log(v[t])+(1/2)*(F[t]/v[t])
    
    Likelihood = np.sum(GARCH_Dens[1:-1])
    return Likelihood

def GARCH_PROD(params, Y0, T):
    mu = params[0]
    omega = np.exp(params[1])
    alpha = (1-np.exp(-params[2])) / (1+np.exp(-params[2]))
    beta = params[3]
    Y = np.zeros(T)
    sigma2 = np.zeros(T)
    Y[0] = Y0
    sigma2[0] = omega/(1-alpha)
    for t in range (1,T):
        sigma2[t] = omega+alpha*((Y[t-1]-mu)**2)+beta*(sigma2[t-1]);
        Y[t] = mu-np.sqrt(sigma2[t])*np.random.normal(0,1,1);
    return Y  

param_star = minimize(GARCH, param0, method = 'BFGS',options = {'xtol':1e-8,'disp':True})

T = 29;
mu = openmean1.mean();
sig  = openmean1.var();
Y=openmean1
param0 = np.array([mu,2.5,0.3,0.5])
param_star = minimize(GARCH, param0, method = 'BFGS',options = {'xtol':1e-8,'disp':True})
Y_GARCH = GARCH_PROD(param_star.x,Y[350],T)
#Y_GARCH = GARCH_PROD(param_star.x,Y[0],T)
timevec = np.linspace(1,T,T)
plt.plot(timevec, Y ,'b', timevec[5:-1], Y_GARCH[5:-1],'r:')

rmse = np.sqrt(np.mean((Y_GARCH-Y)**2))
print(rmse)## 2 week data fitting for Garch and Garch-T

fb=facebook[350:379]
fb.head()
fb.shape

openmean1= (fb['Open']+ fb['Close'])/2
#openmean

def GARCH(Y):
    mu = param0[0]
    omega = np.exp(param0[1])
    alpha = (1-np.exp(-param0[2])) / (1+np.exp(-param0[2]))
    beta = param0[3]
    T = Y.shape[0]
    GARCH_Dens = np.zeros(T)
    sigma2 = np.zeros(T)
    F = np.zeros(T)
    v = np.zeros(T)
    
    for t in range(1,T):
        sigma2[t] = omega+alpha*((Y[t-1]-mu)**2)+beta*(sigma2[t-1]);
        F[t] = Y[t] - mu-np.sqrt(sigma2[t])*np.random.normal(0,1,1);
        v[t] = sigma2[t];
        GARCH_Dens[t] = (1/2)*np.log(2*np.pi)+(1/2)*np.log(v[t])+(1/2)*(F[t]/v[t])
    
    Likelihood = np.sum(GARCH_Dens[1:-1])
    return Likelihood

def GARCH_PROD(params, Y0, T):
    mu = params[0]
    omega = np.exp(params[1])
    alpha = (1-np.exp(-params[2])) / (1+np.exp(-params[2]))
    beta = params[3]
    Y = np.zeros(T)
    sigma2 = np.zeros(T)
    Y[0] = Y0
    sigma2[0] = omega/(1-alpha)
    for t in range (1,T):
        sigma2[t] = omega+alpha*((Y[t-1]-mu)**2)+beta*(sigma2[t-1]);
        Y[t] = mu-np.sqrt(sigma2[t])*np.random.normal(0,1,1);
    return Y  

param_star = minimize(GARCH, param0, method = 'BFGS',options = {'xtol':1e-8,'disp':True})

T = 29;
mu = openmean1.mean();
sig  = openmean1.var();
Y=openmean1
param0 = np.array([mu,2.5,0.3,0.5])
param_star = minimize(GARCH, param0, method = 'BFGS',options = {'xtol':1e-8,'disp':True})
Y_GARCH = GARCH_PROD(param_star.x,Y[350],T)
#Y_GARCH = GARCH_PROD(param_star.x,Y[0],T)
timevec = np.linspace(1,T,T)
plt.plot(timevec, Y ,'b', timevec[5:-1], Y_GARCH[5:-1],'r:')

rmse = np.sqrt(np.mean((Y_GARCH-Y)**2))
print(rmse)
df = facebook[['Open']]

df.reset_index(level=0, inplace=True)

df.columns=['ds','y']<jupyter_output><empty_output><jupyter_text>Simple Moving Average<jupyter_code>rolAvg = df.y.rolling(window=5).mean()

rolAvg2 = df.y.rolling(window=25).mean()

plt.plot(df.ds, df.y, label='FB')

plt.plot(df.ds, rolAvg, label='FB 5 Day', color='orange')

plt.plot(df.ds, rolAvg2, label='FB 25 Day', color='magenta')

plt.legend(loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>Exponential Moving Average<jupyter_code>movAvg = df.y.ewm(span=5, adjust=False).mean()

movAvg2 = df.y.ewm(span=25, adjust=False).mean()

plt.plot(df.ds, df.y, label='FB')

plt.plot(df.ds, movAvg, label='FB 5 Day')

plt.plot(df.ds, movAvg2, label='FB 25 Day')

plt.legend(loc='upper left')

plt.show()<jupyter_output><empty_output><jupyter_text>We can see that using this signal we could have predicted the price trend of FB. When short-term crosses above long-term we get a buy signal. When short-term passes below the longer-term we get a sell signal.  Linear Regression of Simple Moving Average<jupyter_code>movAverage= rolAvg
movAverage.shape
# Sample for X,Y
T = 861

#mu = (facebook['Open '].mean(),facebook['Open '].mean(),facebook['Open '].mean()) 
#cov = [[1,0.75,-0.35],[0.75,1,0.9],[-0.35,0.9,1]]
#F = np.random.multivariate_normal(mu,cov,T)
#Sample for Y,X
X = movAvg
X = np.column_stack([np.ones((T,1)),X])
#X=pd.DataFrame(movAverage[4:])
#print(X)

N = X.shape
#beta = np.array([0.56,2.53,2.05,1.78])
#beta.shape=(N[1],1)
#Y =X@beta+np.random.normal(0,1,(T,1))
Y=facebook['Open'].values
print(Y.shape)
invXX = np.linalg.inv(X.transpose()@X)
beta_hat = invXX@X.transpose()@Y
y_hat = X@beta_hat
residuals = Y-y_hat
sigma2 = (1/T)*residuals.transpose()@residuals

sigma = np.sqrt(sigma2)

#variance - covariance of beta_hat
varcov_beta_hat = (sigma2)*invXX
std_beta_hat = np.sqrt(T*np.diag(varcov_beta_hat))

R_square = 1-(residuals.transpose()@residuals)/(T*np.var(Y))

adj_R_square = 1-(1-R_square)*(T-1)/(T - N[1])

#Testing Coefficents:beta_i
#Null Hypotesis

t_stat = (beta_hat.transpose()-0)/std_beta_hat
p_val = 1-ss.norm.cdf(t_stat)

#Test of joint significance
F_stat= (beta_hat.transpose()@np.linalg.inv(varcov_beta_hat)@beta_hat/N[1])/(residuals.transpose()@residuals/(T-N[1]))
p_val_F= 1 - ss.f.cdf(F_stat,N[1]-1, T-N[1])



plt.plot(facebook['Date'],y_hat,color='magenta')
plt.plot(facebook['Date'],Y, color='black')
plt.show()
rmse = np.sqrt(np.mean((y_hat-Y)**2))
print(rmse)<jupyter_output>2.5522796703913553<jupyter_text>Linear Regression of Exponential Moving Average<jupyter_code>movAverage1= (movAvg)
movAverage1.shape
# Sample for X,Y
T = 861
#mu = (facebook['Open '].mean(),facebook['Open '].mean(),facebook['Open '].mean()) 
#cov = [[1,0.75,-0.35],[0.75,1,0.9],[-0.35,0.9,1]]
#F = np.random.multivariate_normal(mu,cov,T)
#Sample for Y,X
X=movAverage1
X = np.column_stack([np.ones((T,1)),X])

N = X.shape
#beta = np.array([0.56,2.53,2.05,1.78])
#beta.shape=(N[1],1)
#Y =X@beta+np.random.normal(0,1,(T,1))
Y=facebook['Open'].values
print(X)
invXX = np.linalg.inv(X.transpose()@X)
beta_hat = invXX@X.transpose()@Y
y_hat = X@beta_hat
residuals = Y-y_hat
sigma2 = (1/T)*residuals.transpose()@residuals

sigma = np.sqrt(sigma2)

#variance - covariance of beta_hat
varcov_beta_hat = (sigma2)*invXX
std_beta_hat = np.sqrt(T*np.diag(varcov_beta_hat))

R_square = 1-(residuals.transpose()@residuals)/(T*np.var(Y))

adj_R_square = 1-(1-R_square)*(T-1)/(T - N[1])

#Testing Coefficents:beta_i
#Null Hypotesis

t_stat = (beta_hat.transpose()-0)/std_beta_hat
p_val = 1-ss.norm.cdf(t_stat)

#Test of joint significance
F_stat= (beta_hat.transpose()@np.linalg.inv(varcov_beta_hat)@beta_hat/N[1])/(residuals.transpose()@residuals/(T-N[1]))
p_val_F= 1 - ss.f.cdf(F_stat,N[1]-1, T-N[1])
plt.plot(facebook['Date'],y_hat,color='orange')
plt.plot(facebook['Date'],Y, color='black')
plt.show()
rmse = np.sqrt(np.mean((y_hat-Y)**2))
print(rmse)<jupyter_output>2.5522796703913553<jupyter_text>Autoregression<jupyter_code>facebook['y_hat'] = y_hat
facebook.head()
#facebook = facebook.set_index(pd.DatetimeIndex(facebook['Date']))
#facebook.head()
from pandas import Series
from pandas import DataFrame
from pandas import concat
from matplotlib import pyplot
values = DataFrame(facebook['y_hat'].values)
dataframe = concat([values.shift(1), values], axis=1)
dataframe.columns = ['t-1', 't']
result = dataframe.corr()
print(result)
from pandas import Series
from matplotlib import pyplot
from statsmodels.tsa.ar_model import AR


# split dataset
X = facebook['y_hat']
#def AR(X):
train, test = X[1:len(X)-10], X[len(X)-10:]
# train autoregression
model = AR(train)
model_fit = model.fit()
#print('Lag: %s' % model_fit.k_ar)
#print('Coefficients: %s' % model_fit.params)

predictions = model_fit.predict(start=800, end=861, dynamic=False)
#print(len(predictions))
#rmseCalc()
#return predictions
# error = mean_squared_error(test, predictions)
# print('Test MSE: %.3f' % error)
# plot results
#X = facebook['y_hat']
#predictions = AR(X)
pyplot.plot(test)
pyplot.plot(predictions, color='red')
pyplot.show()
model_fit.predict(600,610)<jupyter_output><empty_output><jupyter_text>CAPM<jupyter_code>import pandas_datareader as pdr
from pandas_datareader import data, wb
from datetime import date
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

def capm(start_date, end_date, ticker1, ticker2):
    df = pdr.get_data_yahoo(ticker1, start_date, end_date)
    dfb = pdr.get_data_yahoo(ticker2, start_date, end_date)


    rts = df.resample('M').last()
    rbts = dfb.resample('M').last()
    dfsm = pd.DataFrame({'s_open' : rts['Open'], 'b_open' : rbts['Open']}, index=rts.index)


    dfsm[['s_returns', 'b_returns']] = dfsm[['s_open','b_open']]/dfsm[['s_open','b_open']].shift(1) -1
    dfsm = dfsm.dropna()

    covmat = np.cov(dfsm["s_returns"], dfsm["b_returns"])


    beta = covmat[0,1]/covmat[1,1]
    alpha = np.mean(dfsm["s_returns"])-beta*np.mean(dfsm["b_returns"])


    y = beta * dfsm["b_returns"] + alpha
    SS_res = np.sum(np.power(y - dfsm["s_returns"],2))
    SS_tot = covmat[0,0]*(len(dfsm) - 1) # SS_tot is sample_variance*(n-1)

    r_squared = 1.0 - SS_res/SS_tot
    volatility = np.sqrt(covmat[0,0])
    momentum = np.prod(1+dfsm["s_returns"].tail(12).values) - 1.0


    prd = 12.0 
    alpha = alpha*prd
    volatility = volatility*np.sqrt(prd)

    print ("Beta, alpha, r_squared, volatility, momentum:")
    print (beta, alpha, r_squared, volatility, momentum)

    %matplotlib inline
    fig,ax = plt.subplots(1,figsize=(20,10))
    ax.scatter(dfsm["b_returns"], dfsm['s_returns'], label="Data points")
    beta,alpha = np.polyfit(dfsm["b_returns"], dfsm['s_returns'], deg=1)
    ax.plot(dfsm["b_returns"], beta*dfsm["b_returns"] + alpha, color='red', label="CAPM line")

    
    plt.title('Capital Asset Pricing Model, finding alphas and betas')
    plt.xlabel('Market return $R_m$', fontsize=14)
    plt.ylabel('Stock return $R_i$')

    plt.legend()
    plt.grid(True)
    plt.show()
    return beta, y


beta, y = capm('2016-01-04', '2019-08-30','FB', '^GSPC')<jupyter_output>Beta, alpha, r_squared, volatility, momentum:
1.226229612230391 0.025581471551512983 0.2900380609000528 0.2640902322803233 0.05436074069773733<jupyter_text>SVM<jupyter_code>X = facebook[['Open']]
Y = facebook[['Open']].values
from sklearn.svm import SVR  
svr_model = SVR(kernel='rbf', gamma=0.0005)
result = svr_model.fit(X, Y)
 
y_hat = result.predict(X)
plt.plot(y_hat)
plt.plot(Y)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2045616535994298, 'fraction_numerical': 0.04775481111903065, 'mean_word_length': 4.3423329871928, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '63505', 'license': 'MIT', 'n_tokens_mistral': 7312, 'n_tokens_neox': 6878, 'n_words': 1055, 'path': 'Project/Facebook_Portfolio.ipynb', 'repo_name': 'yangtong951019/INFO7390_Final_Project', 'size': '15433'}