import numpy as np import matplotlib.pyplot as plt from scipy.stats import multivariate_normal # Import stock price data with open('StockPrices_2021.npy', 'rb') as f: t = np.load(f) S = np.load(f) N = S.shape[0] # Amount of data points per realization K = S.shape[1] # Amount of realizations dt = 0.1 # Time step between data points (days)

# Plot the first realization as a function of time. plt.plot( np.arange(0, N*dt, step=dt), S[:,591]) plt.xlabel("time") plt.ylabel("stock price") plt.show() # Plot the first 100 realizations together as a function of time. for i in range(100): plt.plot( np.arange(0, N*dt, step=dt), S[:,i]) plt.xlabel("time") plt.ylabel("stock price") plt.show()

# Plot the mean of the stock prices as a function of time. means = [] for i in range(N): means.append( sum(S[i,:])/K ) plt.plot( np.arange(0, N*dt, step=dt), means) plt.xlabel("time") plt.ylabel("stock price mean") plt.show() # Plot the variance of the stock prices as a function of time. variances = [] print(S[:0,:]) for i in range(N): variances.append(sum((S[i,:]-means[i])**2) / ( K - 1 ) ) plt.plot( np.arange(0, N*dt, step=dt), variances) plt.xlabel("time") plt.ylabel("stock price variance") plt.show()

# Compute the correlation matrix Rs. Rs = np.ones((N,N)) for i in range(0,N): for j in range(0,N): Rs[i,j] = np.mean( S[i,:] * S[j,:] ) # Visualize Rs using plt.matshow(). plt.matshow(Rs) plt.colorbar() plt.xlabel("i") plt.ylabel("j")

# Apply the transformation to the data. mat = [] for j in range(N-1): log_list = np.log( S[j+1,:] / S[j,:] ) mat.append(log_list) mat = np.stack(mat, axis=0) print(mat) #print(mat) #mat = np.asmatrix(mat) # Plot the mean of the transformed realizations as a function of time. #the first index number signifies the timestep 0-100 x_means = [] for i in range(N-1): x_means.append( sum(mat[i,:])/K ) plt.plot( np.arange(0, (N-1)*dt, step=dt), x_means) plt.xlabel("time") plt.ylabel("mean of transformed stock price") plt.show() # Plot the variance of the realizations as a function of time. x_variances = [] for i in range(N-1): #x_variances.append( sum( (mat[i,:]-x_means[i])**2 ) / ( K ) ) x_variances.append(np.var(mat[i,:])) plt.plot( np.arange(0, (N-1)*dt, step=dt), x_variances) plt.xlabel("time") plt.ylabel("variance of transformed stock price") plt.show() # Plot the variance of the realizations as a function of time.

# Compute the correlation matrix Rx. Rx = np.ones((N-1,N-1)) for i in range(0,N-1): for j in range(0,N-1): Rx[i,j] = np.mean( mat[i,:] * mat[j,:] ) # Visualize Rx using plt.matshow(). plt.matshow(Rx) plt.colorbar() plt.xlabel("i") plt.ylabel("j")

#2. Calculate the (constant) mean and variance of the process Xn. mx = sum(x_means)/len(x_means) vx = sum(x_variances)/len(x_variances) # Print the mean and variance. print(f'Mean = {mx:e}') print(f'Variance = {vx:e}')

# Compute the time average for a single realization of Xn. mxt = [] print(mat[:1, 338]) for i in range(N-1): mxt.append( sum( mat[:(i+1),338] ) / (i+1) ) # Plot the time average for a single realization of Xn. plt.plot( np.arange(0, (N-1)*dt, step=dt), mat[:,338], label= "k2=338") plt.plot( np.arange(0, (N-1)*dt, step=dt), mxt, label="running average") plt.xlabel("time") plt.ylabel("transformed stock price ") plt.legend() plt.show() # Print the final value of the time average. print(f'Final value of time average = {mxt[-1]:e}')

# Compute the drift mu and the volatility sigma. sigma = vx / (( dt )**(1/2)) mu = mx / dt + 0.5 * sigma * sigma # Print the drift and volatility. print(f'mu = {mu:e}') print(f'sigma = {sigma:e}')

# Plot the mean of Sn, both from the model and from the data. #sigma=0.15 exp_sn = [] var_sn = [] for n in range(N): exp_sn.append(means[0] * np.exp( mu * n * dt) ) #var_sn.append( means[0] * means[0] * np.exp( 2 * mu * n * dt) * ( np.exp( ( sigma ** 2 )* n * dt ) - 1 )) var_sn.append( means[0]**2 * np.exp( 2 * mu * n * dt) * ( np.exp( sigma*sigma * n * dt) - 1 ) ) plt.plot( np.arange(0, (N)*dt, step=dt), exp_sn, label = "analytic") plt.plot( np.arange(0, (N)*dt, step=dt), means, label = "experimental") plt.legend() plt.xlabel("time") plt.ylabel("stock price mean") plt.show() plt.plot( np.arange(0, N*dt, step=dt), var_sn, label = "analytic") plt.plot( np.arange(0, N*dt, step=dt), variances, label = "experimental") plt.legend() plt.xlabel("time") plt.ylabel("stock price variance") plt.show()

from scipy.stats import multivariate_normal # Use the multivariate_normal object from scipy.stats to fit a multivariate normal distribution on the data. res = 0.001 xy = np.mgrid[-0.15:0.15:res, -0.15:0.15:res].reshape(2,-1).T X1 = mat[91,:] X2 = mat[38,:] fit = multivariate_normal(mean=[sum(X1)/len(X1),sum(X2)/len(X2)], cov = [np.cov(X1),np.cov(X2)]) fit_results = fit.pdf(xy) #print(fit_results) #Plot the 2D histogram of the instances of the transformed data, as well as the multivariate normal distribution fit from scipy plt.hist2d(X1,X2, bins = int(( K )**(1/2)), density=1) plt.colorbar() plt.xlabel("x") plt.ylabel("y") plt.title("2D histogram for 2 time instances") plt.show() plt.scatter(xy[:,0],xy[:,1],c=fit_results) #plt.tricontourf(xy[:,0],xy[:,1], fit_results) plt.xlabel("x") plt.ylabel("y") plt.xlim([-0.14,0.14]) plt.ylim([-0.14,0.14]) plt.title("Fitted normal distribution") plt.colorbar() #plt.scatter(np.arange(-0.15,0.15,0.01), np.arange(-0.15,0.15,0.01),c=fit_results) plt.show() #plt.plot( np.arange(0, (K)*dt, step=dt), X1 ) #plt.plot( np.arange(0, (K)*dt, step=dt), fit_results ) #plt.show()

# For time instance n_a, Use the multivariate_normal object from scipy.stats to fit a 1D normal distribution on the data. sup = np.arange(-0.15,0.15,0.01) # For time instance n_a, plot the histogram of the instances of the transformed data, # as well as the 1D normal distribution fit from scipy. fit = multivariate_normal(mean=sum(X1)/len(X1), cov = [np.cov(X1)], allow_singular=True) fit_results = fit.pdf(np.arange(-0.15,0.15,res)) plt.plot(np.arange(-0.15,0.15,res), fit_results) plt.hist(X1, bins = int(( K )**(1/2)), density = 1) plt.xlabel("value") plt.ylabel("density") plt.title("n_91") plt.show() fit = multivariate_normal(mean=sum(X2)/len(X2), cov = [np.cov(X2)], allow_singular=True) fit_results = fit.pdf(np.arange(-0.15,0.15,res)) plt.plot(np.arange(-0.15,0.15,res), fit_results) plt.hist(X2, bins = int(( K )**(1/2)), density = 1) plt.xlabel("value") plt.ylabel("density") plt.title("n_38") plt.show() # For time instance n_b, Use the multivariate_normal object from scipy.stats to fit a 1D normal distribution on the data. # For time instance n_b, plot the histogram of the instances of the transformed data, # as well as the 1D normal distribution fit from scipy.