#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  7 19:38:47 2019

@author: jmadriga
"""
#imports libraries 
import numpy as np
from scipy.stats import norm

np.random.seed(42) #sets seed for reproducibility

#Matrix
sigma=np.array([[4,-1],[-1,4]])
A=np.linalg.cholesky(sigma)
#number of samples
N=int(1E6)
C_alpha=norm.ppf(1-0.05/2)

# Defines function and the cv function 
psi = lambda x,a: 1.0*np.prod(x>a,1)
psi_cv= lambda x,a: 1.0*(np.sum(x,1)>2*a)

#samples iid from N
X=A@np.random.standard_normal((2,N))

#test values for a

avals=np.arange(1,8)

for i in range(len(avals)):
    #samples iid from N
    X=(A@np.random.standard_normal((2,N))).T #the T is a transpose, so 
    #that it is a matrix of Nx2 samples 
    a=avals[i]
    #creates Z and its cv Y
    Z=psi(X,a)
    Y=psi_cv(X,a)
    #obtaines the sample covariance matrix
    C=np.cov(Z,Y)
    EY=1-norm.cdf(2*a/np.sqrt(6)) #exact mean of Y
    alpha_hat=-C[0,1]/C[1,1]
    rho2=np.corrcoef(Z,Y)[0,1]  #computes the correlation coeff
    mu_cmc=round(np.mean(Z),6)
    Z_tilde=Z+alpha_hat*(Y-EY) #Writes down \tilde Z , which is the correlated Var
    mu_cv=round(np.mean(Z_tilde),6) #computes the mean
    
    #Computes the standar errors of Z and \tilde Z. 
    SE_mc=round(C_alpha*(np.var(Z)/N)**0.5,6)
    SE_cv=round(C_alpha*(np.var(Z_tilde)/N)**0.5,6)
    
    print('----------------------------------------------------')
    print('a= '+str(a))
    print('Crude Monte Carlo mean '+str(mu_cmc)+ ' +- '+str(SE_mc))
    print('control Variate mean   '+str(mu_cv)+  ' +- '+str(SE_cv))
    print('SEmc/SEcv              '+str(SE_mc/SE_cv))