"""
Author: Christopher Ting
Date: 2019-07-10
"""

from __future__  import division, print_function
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt


def indices(a, func):
   return [i for (i, val) in enumerate(a) if func(val)]

def autocorrelation(x):
   """
   http://stackoverflow.com/q/14297012/190597
   http://en.wikipedia.org/wiki/Autocorrelation#Estimation
   """
   n = len(x)
   variance = x.var()
   x = x-x.mean()
   r = np.correlate(x, x, mode = 'full')[-n:]
   assert np.allclose(r, np.array([(x[:n-k]*x[-(n-k):]).sum() for k in range(n)]))
   result = r/(variance*(np.arange(n, 0, -1)))
   return result[1]

#########################################################################################
# Reading in the data file from CRSP

data_dir = '.\\'
ticker = 'GE'

fn = data_dir + ticker + ".csv"
data = pd.read_csv(fn, low_memory=False)

start_date, end_date = 19251231, 20111230

a = data.date.between(start_date, end_date, 'both')
b = indices(a, lambda x: x==True)
data = data.iloc[b]

prc, so = data['PRC'], data['SHROUT']
r, d = data['RETX'],  data['date']
f1, f2 = data['CFACPR'], data['CFACSHR']

n, old_d, c = len(d), 0, 0
tdate, p = ["" for i in range(n)], np.zeros(n, dtype=float)
for i in range(n):
   if d[i] == old_d:
      print(d[i], r[i], r[i-1], prc[i], prc[i-1])
   else:
      if np.isnan(prc[i]) == False and prc[i] > 0:
         old_d = d[i]
         p[c] = prc[i]/f1[i]
         tdate[c] = old_d
         c += 1
      else:
         print(d[i], r[i], r[i-1], prc[i], prc[i-1])

p = p[:c]
tdate = tdate[:c]

log_price = np.log(p)
r1 = np.diff(log_price)

###############################################################################
# Variance ratio test

print('\n')
mr1, sr1 = np.mean(r1), np.var(r1, ddof=0)
VR, Z = np.ones(11, dtype=float), np.zeros(11, dtype=float)
auto = np.zeros(11, dtype=float)

nq, nobs = 11, [0]
T = len(r1)
for q in range(2, nq):
   idx = range(0, c, q)
   p_q = log_price[idx]
   r_q = np.diff(p_q)

   auto[q] = autocorrelation(r_q)

   print('q = %d, length of p_q = %d, length of r_q = %d' % \
         (q, len(p_q), len(r_q)))

   M = len(r_q)
   s = 0
   mrq = q*mr1
   for i in range(M):
      s += (r_q[i] - mrq)**2
   srq = s/M 
   VR[q] = srq/(q*sr1)
   Z[q] = np.sqrt(T)*(VR[q]-1)/np.sqrt(2*(q-1))
   nobs.append(M)

VR, Z, auto = VR[1:q+1], Z[1:q+1],  auto[1:q+1]

q = range(1,11)
plt.plot(q, VR, '-^r')
plt.show()


print('\n')
for i in q[1:]:
   i1 = i-1
   print('%02d %d\t %0.3f\t %0.4f\t %0.2f' % 
         (i, nobs[i1], auto[i1], VR[i1], Z[i1]))