<jupyter_start><jupyter_text>Empirical Analysis * This notebook contains the codes to provide the data and analysis in Section 5.* You can view and download this notebook at my [github repo](https://github.com/okuchap/DAA_Lucas). (https://github.com/okuchap/DAA_Lucas)<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import pathlib
import time
from datetime import datetime
from scipy.optimize import minimize

# used for greedy optimization
from sklearn.utils.extmath import cartesian

#np.set_printoptions(precision=3)
np.set_printoptions(suppress=True)
# plt.style.use("ggplot")
# plt.rcParams["font.size"] = 13
mpl.rcParams["font.family"] = 'Osaka'<jupyter_output><empty_output><jupyter_text>Overview * Focus on the last quarter of 2018 (October 2018 to December 2018)* Assume that the hash supply is linear in reward:$$H(t) = \alpha + \beta R(t)$$The estimates we obtained are$$\hat{\alpha} = 23.60, \hat{\beta} = 9.02, \quad\sigma(\hat{\alpha}) = 2.25, \sigma(\hat{\beta}) = 0.99$$  Method  Data Cleaning * Compute the block time using timestamps.* Basically, the time required to generate the $i$-th block is $\Delta_i := t_i - t_{i-1}$.    - However, sometimes $\Delta_i$ becomes negative due to inaccuracy.* We conduct our analysis as follows:Let $\{t_i\}_i$ be the times recorded in timestamps. For each $i$, let$$k^* = \min_{k \geq 1} \{k \colon t_i < t_{i+k}\}$$and linearly interpolate $t_0, t_1, \dots, t_k$ by the function `adjust_BGT`.<jupyter_code>def adjust_unixtime(raw_times):
    '''
    adjust block generation time.
    Conduct linear interpolation if the order of time stamps are reversed.
    
    Parameters
    ----------
    raw_times:
        pandas Series containing raw time stamps(unixtime)
        The index should starts at zero.
        
    Returns
    -------
    times:
        Adjusted time stamps.
    '''
    times = raw_times.copy()
    for i in range(times.shape[0] - 1):
        if times[i+1] - times[i] < 0:
            count = 1
            diff = -1
            while diff < 0:
                count += 1
                diff = times[i+count] - times[i]
            delta = times[i+count] - times[i]
            for j in range(count-1):
                times[i+j+1] = times[i] + (delta/count)*(j+1)

    return times<jupyter_output><empty_output><jupyter_text>Estimate the hash supply function in BTC system $$h(R) = \alpha + \beta R$$* The loglikelihood is$$\log \mathcal{L}(\alpha,\beta) = \sum_{t} \left\{ \log D(t) + \log(\alpha + \beta R(t)) -D(t)(\alpha + \beta R(t))B(t) \right\}\propto\sum_{t} \left\{\log(\alpha + \beta R(t)) -D(t)(\alpha + \beta R(t))B(t) \right\}$$ * Focus on October 2018 to December 2018. Here, we include the transaction fee in expected reward. * The file `BTCdata_2.csv` contains the raw data.<jupyter_code>path = '../data/BTCdata_2.csv'
df_BTC = pd.read_csv(path, header=0)
df_BTC = df_BTC.rename(columns={'unixtime_precise': 'unixtime'})
df_BTC = df_BTC.loc[: ,['height', 'unixtime', 'probability of success /Eh', 'expected return (USD/Eh)']]
df_BTC['time'] = pd.to_datetime(df_BTC['unixtime'].astype(int), unit='s')
df_BTC = df_BTC.sort_values(by='height', ascending=True)
df_BTC.head()
df_BTC_2018OctDec = df_BTC[(df_BTC['time'] >= datetime(2018,9,30)) & (df_BTC['time'] <= datetime(2019,1,1))]
df_BTC_2018OctDec = df_BTC_2018OctDec.reset_index(drop=True)
BTC_times = df_BTC_2018OctDec['unixtime']
df_BTC_2018OctDec['adj_time'] = adjust_unixtime(BTC_times)
df_BTC_2018OctDec['adj_block_gen_time(min)'] = df_BTC_2018OctDec['adj_time'].diff()/60
df_BTC_2018OctDec = df_BTC_2018OctDec[df_BTC_2018OctDec['time'] >= datetime(2018,10,1)]
df_BTC_2018OctDec.head()
df = df_BTC_2018OctDec.loc[:, ['adj_block_gen_time(min)', 'probability of success /Eh', 'expected return (USD/Eh)']]
df['adj_block_gen_time(sec)'] = df['adj_block_gen_time(min)']*60
df.head()<jupyter_output><empty_output><jupyter_text>Optimization via Scipy<jupyter_code>block_times = df.loc[:, 'adj_block_gen_time(sec)'].values
block_gen_rates = df.loc[:, 'probability of success /Eh'].values
exp_rewards = df.loc[:, 'expected return (USD/Eh)'].values

def log_likelihood(params):
    '''
    params = [a, b]
    '''    
    a = params[0]
    b = params[1]
    
    temp = np.log(a + b * exp_rewards) \
        - (block_gen_rates * (a + b * exp_rewards) * block_times)
    ll = np.sum(temp)
    
    return -ll
x0 = np.array([24.0, 9.0])

result = minimize(fun = log_likelihood,
                  x0=x0, 
                  method="L-BFGS-B")
print(result)<jupyter_output>fun: -35200.9856368449
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([-0.0007276, -0.0007276])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 36
      nit: 10
     njev: 12
   status: 0
  success: True
        x: array([23.58815554,  9.02489233])<jupyter_text>Therefore, the MLE is$$\hat{\alpha} = 23.59, \hat{\beta} = 9.02.$$  Estimation via log-scaled variables * Assume the following linear model and estimate parameters $\alpha, \beta$ via MLE:$$\log H(t) = \alpha + \beta \log R(t),$$which is equivalent to$$H(t) = \mathrm e^{\alpha} R(t)^{\beta}.$$* The likelihood function is$$\mathcal{L} = \prod_t W(t) \left(\mathrm e^{\alpha} R(t)^{\beta}\right)\exp \left(-W(t) B(t) \mathrm e^{\alpha} R(t)^{\beta}\right).$$* The log-likelihood function is$$\log \mathcal{L} \propto \sum_t \left[\alpha + \beta \log R(t)  - W(t) B(t) \mathrm e^{\alpha} R(t)^{\beta} \right].$$ We obtained the MLE such that$$\hat \alpha = 3.42, \ \hat \beta = 0.45.$$  Computation<jupyter_code>def log_likelihood_2(params):
    '''
    params = [a, b]
    '''    
    a = params[0]
    b = params[1]
    
    temp = a + b*np.log(exp_rewards) - block_gen_rates*block_times*np.exp(a)*np.power(exp_rewards, b)
    ll = np.sum(temp)
    
    return -ll
x0 = np.array([3.0, 0.5])

result = minimize(fun = log_likelihood_2,
                  x0=x0, 
                  bounds = ((0,np.inf),(0.01,np.inf)),
                  method="L-BFGS-B")
print(result)<jupyter_output>fun: -35200.493577691355
 hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>
      jac: array([0.00800355, 0.00582077])
  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
     nfev: 39
      nit: 11
     njev: 13
   status: 0
  success: True
        x: array([3.41878738, 0.45145085])<jupyter_text>Therefore, the MLE is$$\hat \alpha = 3.42, \ \hat \beta = 0.45.$$  Plot<jupyter_code>@np.vectorize
def log_hash_supply(log_x):
    return np.log(55) - np.log(1 + np.exp(4.5)*np.exp(-3*np.exp(log_x)))
x = np.linspace(0.01, 3.0, 300)
log_x = np.log(x)
alpha = 3.42
beta = 0.45
y = log_hash_supply(log_x)
y_est = alpha + beta * log_x
flag = log_x > np.log(1.58)
log_x_truncated = log_x[flag]
y_est_truncated = y_est[flag]
plt.rcParams["font.size"] = 12

fig = plt.figure(facecolor='w')
ax = fig.add_subplot(1,1,1, aspect='equal', title='Hash Supply Function')
ax.plot(log_x, y, label='Assumed')
ax.plot(log_x_truncated, y_est_truncated, label='Estimated', linestyle='--')
ax.vlines(np.log(1.58), ymin=-0.5, ymax=4.05, linestyles='-', colors='gray')

ax.set_xlabel('Log Reward Rate (USD/Ehash)')
ax.set_ylabel('Log Hash Rate (EHash/second)')

ax.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Variance-Covariance Matrix The Fisher information matrix is$$I(\hat \alpha, \hat \beta) = \begin{bmatrix}\sum_t 1 & \sum_t \log R(t) \\\sum_t \log R(t) & \sum_t (\log R(t))^2\end{bmatrix}.$$Thus the variance-covariance matrix is$$V(\hat \alpha, \hat \beta) = (I (\hat \alpha, \hat \beta))^{-1}=\begin{bmatrix}0.00223838 & -0.00263739 \\-0.00263739 & 0.00322153\end{bmatrix}.$$<jupyter_code>I = np.zeros((2,2))
df.head()
logR = np.log(df['expected return (USD/Eh)'].values)
logR
logR**2
T = df.shape[0]
I[0,0] = T
I[0,1] = np.sum(logR)
I[1,0] = np.sum(logR)
I[1,1] = np.sum(logR**2)
I
V = np.linalg.inv(I)
V<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16924820551567812, 'fraction_numerical': 0.048356630147336606, 'mean_word_length': 3.7728365384615383, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 32, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '525819', 'license': 'MIT', 'n_tokens_mistral': 3347, 'n_tokens_neox': 3047, 'n_words': 679, 'path': 'notebook/empirical_analysis.ipynb', 'repo_name': 'okuchap/DAA_Lucas', 'size': '7941'}