In [216]:
from IPython import parallel
clients = parallel.Client(profile='parallel')
In [222]:
print clients.ids
print "Total %i cores"%(len(clients.ids))
In [223]:
%%px --local
import sys
sys.path.append("\\\\DAP-NAS\\work\\CNOP")
import cPickle as pickle
import numpy as np
import pandas as pd
from statsmodels.tsa.arima_process import arma_generate_sample
import matplotlib.pyplot as plt
import numpy as np
import CNOP
CNOP = reload(CNOP)
CNOP = CNOP.CNOP
import winsound
def threshold(x,thresholds=[],values=[-1,0,1]):
for threshold,val in zip(thresholds,values):
if x < threshold:
return val
return values[-1]
import time
from itertools import repeat
import os
from datetime import datetime
import numpy as np
import numpy.linalg as linalg
from statsmodels.tsa.arima_process import arma_generate_sample
import warnings
In [224]:
%%px --local
x2000 = pd.read_csv("\\\\DAP-NAS\\work\\CNOP\\x2000.csv", delimiter=";")
In [225]:
%%px --local
N=100
def autocorrelated_errors(rho, df, give_distortion=True):
###########################################################
##### This code is for partial-overlap case of CNOP########
###########################################################
N=500
#PreSampling:
df = df[:N]
#df = x2000[:N]
### DGP generation
regime = pd.DataFrame()
beta, alpha = [0.6, 0.4], [0.9, 1.5]
gammam, mum = [0.3, 0.9], [-0.67, 0.36]
gammap, mup = [0.2, 0.3], [0.02, 1.28]
distortion = np.array([arma_generate_sample(ar = [1, -rho], ma=[1], nsample = N) for i in range(3)])
regime["xbeta"] = df[["X1", "X2"]].dot(beta) + distortion[0]
regime['z-gammam'] = df[["X1", "X3"]].dot(gammam) + distortion[1]
regime['z+gammap'] = df[["X2", "X3"]].dot(gammap) + distortion[2]
regime['regime'] = regime['xbeta'].apply(lambda x: threshold(x,thresholds=alpha))
regime['Y-']=regime['z-gammam'].apply(lambda x: threshold(x, thresholds=mum,values=[-2,-1,0]))
regime['Y+']=regime['z+gammap'].apply(lambda x: threshold(x, thresholds=mup,values=[0,1,2]))
df['Y'] = 0
df['Y'] += np.where(regime['regime']==-1,regime['Y-'],0)
df['Y'] += np.where(regime['regime']==1,regime['Y+'],0)
###df is full data matrix
#Model starts here:
exog = df[["X1", "X2", "X3"]]
endog = df[["Y"]]
l = {0:df}
pan = pd.Panel(l)
y = pan.ix[:,:,['Y']]
x = pan.ix[:,:,["X1", "X2"]]
zminus = pan.ix[:,:,["X1", "X3"]]
zplus = pan.ix[:,:,["X2", "X3"]]
exog = {'x':x,'zplus':zplus,'zminus':zminus}
CNOP4 = CNOP(y,exog, model='CNOP',interest_step=1)
#Counting execution time for fit ...
exec_time = {"fit":0,"se":0}
start_time = time.time()
res = CNOP4.fit( maxiter=250, disp=0)
exec_time["fit"] = time.time() - start_time
# And for Standard Errors as well
try:
start_time = time.time()
res["se"] = CNOP4.se(res.x)
exec_time["se"] = time.time() - start_time
res['status'] = "OK"
except Exception, e:
print e
res['status'] = e
res["exec_time"]=exec_time
if give_distortion:
return res, distortion
else:
return res
In [123]:
%matplotlib inline
In [146]:
rho=0.5
res, distortion = autocorrelated_errors(rho=rho, df=x2000)
print "Done for rho=%.3f" % (rho)
print "Res.status:", res.success
beta, alpha = [0.6, 0.4], [0.9, 1.5]
gammam, mum = [0.3, 0.9], [-0.67, 0.36]
gammap, mup = [0.2, 0.3], [0.02, 1.28]
res_real = beta+alpha+gammam+mum+gammap+mup
plt.bar(np.arange(len(res.x-res_real)),
list(res.x-res_real));
winsound.PlaySound(u'C:\Windows\Media\Windows Print complete.wav', winsound.SND_FILENAME)
print "X & SE:"
print np.dstack((res.x,res.se))
In [226]:
%%px --local
def runMC(rho, n_repl, fun, append_folder=None, path="\\\\DAP-NAS\\work\\CNOP\\dumps\\", prefix = "CNOPres", ext = ".p", **kwargs):
print "Starting MC on %i cores, %i replications of \"%s\" with rho = %.2f observations"\
%(len(clients.ids), n_repl, fun.__name__, rho)
print "If you again run next(), you will get wait_interactive() form and results would be backed up"
print "Untill then you are free to use AsyncResult object that was apready yielded, but data will not be backed up!"
view = clients.load_balanced_view()
def doJob(i):
rho,path,fun, kwargs=i
dumb_str=time.strftime("%Y%m%d%H%M%S")+str(np.random.rand())
filename = path + dumb_str
try:
local_res = fun(rho, **kwargs)
pickle.dump(local_res, open( filename, "wb" ) )
return local_res, filename
except Exception, e:
pickle.dump(e, open( filename, "wb" ) )
return Exception, e, filename
if append_folder is None:
cur_time_str=time.strftime('%Y%m%d_%H%M%S')
temp_folder_str = path+'temp/'+cur_time_str+"/"
if not os.path.exists(temp_folder_str):
os.makedirs(temp_folder_str)
else:
temp_folder_str = path+'temp/'+append_folder+"/"
print temp_folder_str
if append_folder is not None:
readme_f = file(temp_folder_str+"!README.txt", "w")
readme_f.write("Doing MC on %i cores, %i replications of \"%s\" with rho = %.2f observations"\
%(len(clients.ids), n_repl, fun.__name__, rho))
readme_f.close()
ar = view.map_async(doJob, [[rho, temp_folder_str, fun, kwargs]]*n_repl)
yield ar
ar.wait_interactive()
results = ar.get()
cur_time_str=time.strftime('%Y%m%d_%H%M%S')
filename = path + prefix + cur_time_str + ext
print "DONE! DB in %s"%(temp_folder_str)
yield results
In [66]:
MCgenr03 = runMC(.3, 10000, autocorrelated_errors, df=x2000)
MCAsyncMapResult03 = next(MCgenr03)
In [67]:
results03=next(MCgenr03)
In [68]:
MCAsyncMapResult03.abort()
In [74]:
i=MCAsyncMapResult03.get()
In [80]:
someresults = MCAsyncMapResult03
In [ ]:
In [86]:
MCgenr03 = runMC(.3, 9300, autocorrelated_errors, df=x2000)
MCAsyncMapResult03 = next(MCgenr03)
In [87]:
results03=next(MCgenr03)
In [89]:
MCAsyncMapResult03.abort()
In [90]:
i=MCAsyncMapResult03.get()
In [92]:
someresults = MCAsyncMapResult03
In [ ]:
In [100]:
MCgenr03 = runMC(.3, 9000, autocorrelated_errors, df=x2000)
MCAsyncMapResult03 = next(MCgenr03)
In [101]:
results03=next(MCgenr03)
In [ ]:
In [112]:
MCgenr06 = runMC(.6, 10000, autocorrelated_errors, df=x2000)
MCAsyncMapResult06 = next(MCgenr06)
In [113]:
results06=next(MCgenr06)
In [ ]:
In [114]:
MCgenr09 = runMC(.9, 10000, autocorrelated_errors, df=x2000)
MCAsyncMapResult09 = next(MCgenr09)
In [115]:
results09=next(MCgenr09)
In [116]:
MCAsyncMapResult09.abort()
In [ ]:
In [129]:
MCgenr09_1 = runMC(.9, 5163, append_folder="20150413_112135", fun=autocorrelated_errors, df=x2000)
MCAsyncMapResult09_1 = next(MCgenr09_1)
In [130]:
results09_1=next(MCgenr09_1)
In [227]:
%%px --local
def autocorrelated_errors_x(rho, df, give_distortion=True):
###########################################################
##### This code is for partial-overlap case of CNOP########
###########################################################
N=500
#PreSampling:
df = df[:N]
#df = x2000[:N]
### DGP generation
regime = pd.DataFrame()
beta, alpha = [0.6, 0.4], [0.9, 1.5]
gammam, mum = [0.3, 0.9], [-0.67, 0.36]
gammap, mup = [0.2, 0.3], [0.02, 1.28]
distortion = arma_generate_sample(ar = [1, -rho], ma=[1], nsample = N)
distortion = np.vstack((distortion, np.random.randn(2,N)))
regime["xbeta"] = df[["X1", "X2"]].dot(beta) + distortion[0]
regime['z-gammam'] = df[["X1", "X3"]].dot(gammam) + distortion[1]
regime['z+gammap'] = df[["X2", "X3"]].dot(gammap) + distortion[2]
regime['regime'] = regime['xbeta'].apply(lambda x: threshold(x,thresholds=alpha))
regime['Y-']=regime['z-gammam'].apply(lambda x: threshold(x, thresholds=mum,values=[-2,-1,0]))
regime['Y+']=regime['z+gammap'].apply(lambda x: threshold(x, thresholds=mup,values=[0,1,2]))
df['Y'] = 0
df['Y'] += np.where(regime['regime']==-1,regime['Y-'],0)
df['Y'] += np.where(regime['regime']==1,regime['Y+'],0)
###df is full data matrix
#Model starts here:
exog = df[["X1", "X2", "X3"]]
endog = df[["Y"]]
l = {0:df}
pan = pd.Panel(l)
y = pan.ix[:,:,['Y']]
x = pan.ix[:,:,["X1", "X2"]]
zminus = pan.ix[:,:,["X1", "X3"]]
zplus = pan.ix[:,:,["X2", "X3"]]
exog = {'x':x,'zplus':zplus,'zminus':zminus}
CNOP4 = CNOP(y,exog, model='CNOP',interest_step=1)
#Counting execution time for fit ...
exec_time = {"fit":0,"se":0}
start_time = time.time()
res = CNOP4.fit( maxiter=250, disp=0)
exec_time["fit"] = time.time() - start_time
# And for Standard Errors as well
try:
start_time = time.time()
res["se"] = CNOP4.se(res.x)
exec_time["se"] = time.time() - start_time
res['status'] = "OK"
except Exception, e:
print e
res['status'] = e
res["exec_time"]=exec_time
if give_distortion:
return res, distortion
else:
return res
In [147]:
rho=0.5
res, distortion = autocorrelated_errors_x(rho=rho, df=x2000)
print "Done for rho=%.3f" % (rho)
print "Res.status:", res.success
beta, alpha = [0.6, 0.4], [0.9, 1.5]
gammam, mum = [0.3, 0.9], [-0.67, 0.36]
gammap, mup = [0.2, 0.3], [0.02, 1.28]
res_real = beta+alpha+gammam+mum+gammap+mup
plt.bar(np.arange(len(res.x-res_real)),
list(res.x-res_real));
winsound.PlaySound(u'C:\Windows\Media\Windows Print complete.wav', winsound.SND_FILENAME)
print "X & SE:"
print np.dstack((res.x,np.array(res.se)))
In [212]:
MCgenr03 = runMC(.3, 8503, fun=autocorrelated_errors_x, df=x2000)
MCAsyncMapResult03 = next(MCgenr03)
In [213]:
results03=next(MCgenr03)
In [ ]:
In [220]:
MCgenr06 = runMC(.6, 10000, fun=autocorrelated_errors_x, df=x2000)
MCAsyncMapResult06 = next(MCgenr06)
In [221]:
results06=next(MCgenr06)
In [ ]:
In [231]:
results09=next(MCgenr09)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [65]:
(np.matrix([[1,2,3],[4,5,6],[7,8,9],[10,12,12]])-[1,2,3]).mean(axis=1)
Out[65]:
In [17]:
rho=0.5
N=1000
In [18]:
ACrng = np.array([arma_generate_sample(ar = [1, -rho], ma=[1], nsample = N) for i in range(3)])
In [33]:
ACrng[0].mean()
Out[33]:
In [36]:
import statsmodels.tsa.stattools
In [39]:
statsmodels.tsa.stattools.acf(ACrng[0])
Out[39]: