Simulating data from the model¶
Description: Random variables induced by the model. Coin flipping example. Simulating consumption and wealth paths from the consumption-savings model.
Stochastic model as random variable¶
$$ \tilde{X} \sim F_{\tilde{X}}(x) $$- distributions directly available in modules
random,NumPy.randomandSciPy.stat - inverse transform sampling to simulated from any cdf
- $ M(\cdot) $ is the structure of the model
- $ \tilde{\varepsilon} $ is the vector of random components in the model (“shocks”)
Example: simulating Markov chains¶
- $ M(\cdot) $ is given by the transition probability matrix $ P $
- $ \tilde{\varepsilon} $ can be though of the vector of standard uniform random variables $ U[0,1] $
- for given initial state $ x_0 $, use the appropriate row of the matrix $ P $ to simulate $ x_1 $
- for given state $ x_i $, use the appropriate row of the matrix $ P $ to simulate $ x_{i+1} $, and continue so forth for a given number of steps
(see video 20)
More interesting example: coin flipping game¶
- fair coin is flipped sequentially
- if 3 heads come up in a row, pay $10
- if 3 tails come up in a row, get $10
What is the distribution of the wins after 100 coin flips?
What is the expected number of tosses to win or loose $100?
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = [12, 8]
In [2]:
def game(p=0.5,nsteps=10,stopping=None,verbose=False):
'''3 tails coin flipping game'''
i, tails, heads, balance = 1, 0, 0, [0.] # initialize all indicators and counters
while True:
coin = np.random.random() > p # draw from U[0,1)
if verbose:
print('T' if coin else 'H',end=' ')
# increment the counters of consecutive heads/tails
tails = tails + 1 if coin else 0
heads = heads + 1 if not coin else 0
i+=1
if tails == 3 or heads == 3:
increment = 10 if tails == 3 else -10
balance.append(balance[-1] + increment)
heads, tails = 0, 0 # four heads or tails should not count two series of three
if verbose:
print('($%1.0f)'%balance[-1],end=' ')
else:
balance.append(balance[-1])
if nsteps and i >= nsteps:
if verbose:
print('Done after %d steps'%nsteps)
break
if stopping and (balance[-1] <= stopping[0] or balance[-1] >= stopping[1]):
if verbose:
print('Done after hitting the boundary')
break
return balance[-1], balance, i
In [3]:
b,bb,st = game(verbose=True)
# b,bb,st = game(verbose=True,stopping=(-30,100),nsteps=None)
print('Balance $%1.1f after %d coin flips'%(b,st))
T T T ($10) T H T T T ($20) H Done after 10 steps Balance $20.0 after 10 coin flips
In [4]:
data = [game(nsteps=250) for i in range(10)]
fig, ax = plt.subplots(figsize=(10,6))
for balance,trajectory,rounds in data:
plt.plot(range(rounds),trajectory,alpha=0.75)
plt.show()
In [5]:
data = [game(nsteps=100) for i in range(5000)]
balances = [i for i,j,k in data]
plt.hist(balances,bins=25,density=True,color='skyblue',edgecolor='k')
plt.title('Distribution of wins/losses after 100 steps of the coin flipping game')
plt.show()
In [6]:
data = [game(nsteps=None,stopping=(-100,200)) for i in range(1000)]
rounds = [k for i,j,k in data]
plt.hist(rounds,bins=25,density=True,color='skyblue',edgecolor='k')
plt.title('Distribution of the number of steps to hit one of the (-100,200) bounds in the coin flipping game')
plt.show()
Simulating from dynamic models¶
General form of the Bellman equation/operator
$$ V(\text{state}) = \max_{\text{decisions}} \big[ U(\text{state},\text{decision}) + \beta \mathbb{E}\big\{ V(\text{next state}) \big| \text{state},\text{decision} \big\} \big] $$Without loss of generality express the next period state is a function of shocks $ g(\text{state},\text{decision},\text{shocks}) = g(x,d,\tilde{\varepsilon}) $
$$ V(x) = \max_{d} \big[ U(x,d) + \beta \mathbb{E}_{\varepsilon}\big\{ V\big( g(x,d,\tilde{\varepsilon}) \big) \big\} \big] = T(V)(x) $$Solution of the model and simulation from the model¶
- Solution for the dynamic model: value function $ V^\star(x) $ satisfying the fixed point condition $ V(x)=T(V)(x) $ AND the policy function
- After the model is solved for $ d^\star(x) $, $ \tilde{X} = M(\tilde{\varepsilon}), \tilde{\varepsilon} \sim F_{\varepsilon}(x) $ applies
- $ \tilde{X} = \{x_0,\tilde{x}_1,\dots,\tilde{x}_T\} $ is the simulated state trajectory (for a given $ x_0 $ and $ T $)
- $ \tilde{\varepsilon} = \{\tilde{\varepsilon}_1,\dots,\tilde{\varepsilon}_T\} $ is a series of shocks, and it holds
The difference between solution and simulation from the dynamic model¶
- Solution: fixed point of the Bellman operator value function $ V^\star(x) $ and the policy function $ d^\star(x) $
- Simulation: a sequence of $ T $ realized points in the state space starting from the initial point $ x_0 $ under the optimal policy
Simulating data from the stochastic consumption-savings model¶
Write the simulator (simulation function) for the stochastic consumption-savings model for:
- given vector of $ N $ initial state points $ M_0 = (M^1_0,\dots,M^N_0) $
- simulate optimal consumption choices and state transitions on each of the $ N $ trajectories for $ T $ periods
- make simulations repeatable by fixing the seed of the random number generator
Simulator algorithm¶
- Solve the model to find the policy function
- Initialize the simulated paths, fix the seed if required
- Loop for $ T $ periods, computing the optimal consumption choices and simulating random incomes
In [7]:
# Developed in the video
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from scipy.stats import lognorm
class deaton():
'''Implementation of the stochastic Deaton consumption-savings problem with random income.'''
def __init__(self, Mbar=10,
ngrid=50, nchgrid=100, nquad=10, interpolation='linear',
beta=.9, R=1.05, sigma=1.):
'''Object creator for the stochastic consumption-savings model'''
self.beta = beta # Discount factor
self.R = R # Gross interest
self.sigma = sigma # Param in log-normal income distribution
self.Mbar = Mbar # Upper bound on wealth
self.ngrid = ngrid # Number of grid points in the state space
self.nchgrid = nchgrid # Number of grid points in the decision space
self.nquad = nquad # Number of quadrature points
self.interpolation = interpolation # type of interpolation, see below
# state and choice space grids, as well as quadrature points and weights are set with setter functions below
def __repr__(self):
'''String representation for the model'''
return 'Deaton model with beta={:1.3f}, sigma={:1.3f}, gross return={:1.3f}\nGrids: state {} points up to {:1.1f}, choice {} points, quadrature {} points\nInterpolation: {}\nThe model is {}solved.'\
.format(self.beta,self.sigma,self.R,self.ngrid,self.Mbar,self.nchgrid,self.nquad,self.interpolation,'' if hasattr(self,'solution') else 'not ')
@property
def ngrid(self):
'''Property getter for the ngrid parameter'''
return self.__ngrid
@ngrid.setter
def ngrid(self,ngrid):
'''Property setter for the ngrid parameter'''
self.__ngrid = ngrid
epsilon = np.finfo(float).eps # smallest positive float number difference
self.grid = np.linspace(epsilon,self.Mbar,ngrid) # grid for state space
@property
def nchgrid(self):
'''Property getter for the nchgrid parameter'''
return self.__nchgrid
@nchgrid.setter
def nchgrid(self,nchgrid):
'''Property setter for the nchgrid parameter'''
self.__nchgrid = nchgrid
epsilon = np.finfo(float).eps # smallest positive float number difference
self.chgrid = np.linspace(epsilon,self.Mbar,nchgrid) # grid for state space
@property
def sigma(self):
'''Property getter for the sigma parameter'''
return self.__sigma
@sigma.setter
def sigma(self,sigma):
'''Property setter for the sigma parameter'''
self.__sigma = sigma
self.__quadrature_setup() # update quadrature points and weights
@property
def nquad(self):
'''Property getter for the number of quadrature points'''
return self.__nquad
@nquad.setter
def nquad(self,nquad):
'''Property setter for the number of quadrature points'''
self.__nquad = nquad
self.__quadrature_setup() # update quadrature points and weights
def __quadrature_setup(self):
'''Internal function to set up quadrature points and weights,
depends on sigma and nquad, therefore called from the property setters
'''
try:
# quadrature points and weights for log-normal distribution
self.quadp,self.quadw = np.polynomial.legendre.leggauss(self.__nquad) # Gauss-Legendre for [-1,1]
self.quadp = (self.quadp+1)/2 # rescale to [0,1]
self.quadp = lognorm.ppf(self.quadp,self.__sigma) # inverse cdf
self.quadw /= 2 # rescale weights as well
self.quadp.shape = (1,1,self.__nquad) # quadrature points in third dimension of 3-dim array
except(AttributeError):
# when __nquad or __sigma are not yet set
pass
def utility(self,c):
'''Utility function'''
return np.log(c)
def next_period_wealth(self,M,c,y=None):
'''Next period budget, returned for all quadrature points'''
M1 = self.R*(M-c) # next period wealth without income
if np.any(y==None):
M1.shape += (1,)*(3-M1.ndim) # add singular dimensions up to 3
# interpolating over income ==> replace with quadrature points
if self.nquad>1:
return M1 + self.quadp # 3-dim array
else:
return M1 # special case of no income
else: # needed for simulated income
if self.nquad>1:
return M1 + y
else:
return M1 # special case of no income
def interp_func(self,x,f):
'''Returns the interpolation function for given data'''
if self.interpolation=='linear':
return interpolate.interp1d(x,f,kind='slinear',fill_value="extrapolate")
elif self.interpolation=='quadratic':
return interpolate.interp1d(x,f,kind='quadratic',fill_value="extrapolate")
elif self.interpolation=='cubic':
return interpolate.interp1d(x,f,kind='cubic',fill_value="extrapolate")
elif self.interpolation=='polynomial':
p = np.polynomial.polynomial.polyfit(x,f,self.ngrid_state-1)
return lambda x: np.polynomial.polynomial.polyval(x,p)
else:
print('Unknown interpolation type')
return None
def bellman_discretized(self,V0):
'''Bellman operator with discretized choice,
V0 is 1-dim vector of values on the state grid
'''
states = self.grid[np.newaxis,:] # state grid as row vector
choices = self.chgrid[:,np.newaxis] # choice grid as column vector
M = np.repeat(states,self.nchgrid,axis=0) # current wealth, state grid in columns (repeated rows)
c = np.repeat(choices,self.ngrid,axis=1) # choice grid in rows (repeated columns)
c *= self.grid/self.Mbar # scale values of choices to ensure c<=M
M1 = self.next_period_wealth(M,c) # 3-dim array with quad point in last dimension
inter = self.interp_func(self.grid,V0) # interpolating function for next period value function
V1 = inter(M1) # value function at next period wealth, 3-dim array
EV = np.dot(V1,self.quadw) # expected value function, 2-dim matrix
MX = self.utility(c) + self.beta*EV # maximand of Bellman equation, 2-dim matrix
MX[c>M] = -np.inf # infeasible choices should have -inf (just in case)
V1 = np.amax(MX,axis=0,keepdims=False) # optimal choice as maximum in every column, 1-dim vector
c1 = c[np.argmax(MX,axis=0),range(self.ngrid)] # choose the max attaining levels of c
return V1, c1
def solve_vfi (self,maxiter=500,tol=1e-4,callback=None):
'''Solves the model using value function iterations (successive approximations of Bellman operator)
Callback function is invoked at each iteration with keyword arguments.
'''
V0 = self.utility(self.grid) # on first iteration assume consuming everything
for iter in range(maxiter):
V1,c1 = self.bellman_discretized(V0)
err = np.amax(np.abs(V1-V0))
if callback: callback(iter=iter,model=self,value=V1,policy=c1,err=err) # callback for making plots
if err < tol:
break # converged!
V0 = V1 # prepare for the next iteration
else: # when iter went up to maxiter
raise RuntimeError('No convergence: maximum number of iterations achieved!')
self.solution = {'value':V1,'policy':c1} # save the model solution to the object
return V1,c1
def solve_plot(self,solver,**kvarg):
'''Illustrate solution
Inputs: solver (string), and any inputs to the solver
'''
if solver=='vfi':
solver_func = self.solve_vfi
else:
raise ValueError('Unknown solver label')
fig1, (ax1,ax2) = plt.subplots(1,2,figsize=(14,8))
ax1.grid(b=True, which='both', color='0.65', linestyle='-')
ax2.grid(b=True, which='both', color='0.65', linestyle='-')
ax1.set_title('Value function convergence with %s'%solver)
ax2.set_title('Policy function convergence with %s'%solver)
ax1.set_xlabel('Wealth, M')
ax2.set_xlabel('Wealth, M')
ax1.set_ylabel('Value function')
ax2.set_ylabel('Policy function')
def callback(**kwargs):
print('|',end='')
grid = kwargs['model'].grid
v = kwargs['value']
c = kwargs['policy']
ax1.plot(grid[1:],v[1:],color='k',alpha=0.25)
ax2.plot(grid,c,color='k',alpha=0.25)
V,c = solver_func(callback=callback,**kvarg)
# add solutions
ax1.plot(self.grid[1:],V[1:],color='r',linewidth=2.5)
ax2.plot(self.grid,c,color='r',linewidth=2.5)
plt.show()
return V,c
def simulator(self,init_wealth=1,T=10,seed=None,plot=True):
'''Simulation of the model for given number of periods from given initial conditions'''
assert hasattr(self,'solution'), 'Need to solve the model before simulating!'
if seed!=None:
np.random.seed(seed) # fix the seed if needed
init_wealth = np.asarray(init_wealth).ravel() # flat np array of initial wealth
N = init_wealth.size # number of trajectories to simulate
sim = {'M':np.empty((N,T+1)),'c':np.empty((N,T+1))}
sim['M'][:,0] = init_wealth # initial wealth in the first column
inter = self.interp_func(self.grid,self.solution['policy']) # interpolation function for policy function
for t in range(T+1):
sim['c'][:,t] = inter(sim['M'][:,t]) # optimal consumption in period t
if t<T:
y = lognorm.rvs(self.sigma,size=N) # draw random income
sim['M'][:,t+1] = self.next_period_wealth(sim['M'][:,t],sim['c'][:,t],y) # next period wealth
if plot:
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(12,6))
ax1.set_title('Simulated wealth and consumption trajectories')
ax1.set_ylabel('Wealth')
ax2.set_ylabel('Consumption')
ax2.set_xlabel('Time period in the simulation')
for ax in (ax1,ax2):
ax.grid(b=True, which='both', color='0.95', linestyle='-')
for i in range(N):
ax1.plot(sim['M'][i,:],alpha=0.75)
ax2.plot(sim['c'][i,:],alpha=0.75)
plt.show()
return sim # return simulated data
In [8]:
m = deaton(ngrid=100,nchgrid=250,sigma=.2,R=1.05,beta=.8,nquad=10)
print(m)
v,c = m.solve_plot(solver='vfi')
sims = m.simulator(init_wealth=[1,2,3],T=5)
Deaton model with beta=0.800, sigma=0.200, gross return=1.050 Grids: state 100 points up to 10.0, choice 250 points, quadrature 10 points Interpolation: linear The model is not solved. ||||||||||
Changes in the class deaton()¶
- save the solution in the object, update the
__repr__()method to indicate if the model is solved - make sure the next period budget can be computed with simulated income
- new
simulator()method - plotting the simulations
In [9]:
# Previously developed simulator code (that requires the changes to the deaton() class!)
def simulator(self, init_wealth=1, T=10, seed=None, plot=True):
'''Simulator for T periods of the model, init_wealth = initial wealth'''
if self.solution is None:
raise(RuntimeError('Model has to be solved before simulations'))
if seed:
np.random.seed(seed) # fix the seed in needed
interp = self.interp_func(self.grid,self.solution['policy']) # interpolator of the policy function
init_wealth = np.asarray(init_wealth).ravel() # initial wealth as flat np array
N = init_wealth.size # number of individuals to simulate
sim = {'M': np.empty((N,T+1),dtype='float'),
'c': np.empty((N,T+1),dtype='float'),
'y': np.empty((N,T+1),dtype='float'),
} # dictionary to hold the simulated values
sim['M'][:,0] = init_wealth # initialize wealth in first period
sim['y'][:,0] = np.full(N,np.NaN) # no income in the first period
for t in range(T+1):
sim['c'][:,t] = interp(sim['M'][:,t]) # optim consumption
if t < T: # until the last period
sim['y'][:,t+1] = np.random.lognormal(sigma=self.sigma,size=N)
sim['M'][:,t+1] = self.next_period_wealth(sim['M'][:,t],sim['c'][:,t],sim['y'][:,t+1])
self.sims = sim # save simulation into the object
if plot:
N = self.sims['M'].shape[0] # number of simulated agents
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(12,6))
ax1.set_title('Simulated wealth and consumption trajectories')
ax1.set_ylabel('Wealth')
ax2.set_ylabel('Consumption')
ax2.set_xlabel('Time period in the simulation')
for ax in (ax1,ax2):
ax.grid(b=True, which='both', color='0.95', linestyle='-')
for i in range(N):
ax1.plot(sim['M'][i,:],alpha=0.75)
ax2.plot(sim['c'][i,:],alpha=0.75)
plt.show()
return sim
m = deaton(ngrid=100,nchgrid=250,sigma=.2,R=1.05,beta=.85,nquad=10)
v,c = m.solve_vfi()
print(m)
# generate initial wealth
init_wealth = np.random.lognormal(sigma=.1,size=10)
# init_wealth = .85*np.ones(shape=5)
simulator(m,init_wealth=init_wealth,T=50);
# print('\nSimulated wealth:\n',m.sims['M'])
# print('\nSimulated consumption:\n',m.sims['c'])
# print('\nIncome realizations:\n',m.sims['y'])
Deaton model with beta=0.850, sigma=0.200, gross return=1.050 Grids: state 100 points up to 10.0, choice 250 points, quadrature 10 points Interpolation: linear The model is solved.
More simulations¶
- Longer horizon, many people
- Change sigma parameter in the model
- Observe tendency to consumption smoothing
In [10]:
m = deaton(R=1.05,beta=.925,nquad=7,sigma=0.5)
m.solve_vfi(maxiter=500)
# m.simulator(T=100,init_wealth=np.random.lognormal(sigma=0.1,size=3))
simulator(m,T=100,init_wealth=np.random.lognormal(sigma=0.1,size=3),plot=True)
Out[10]:
{'M': array([[0.96416167, 1.21328913, 1.02602714, 1.45947934, 1.32219227,
1.69577937, 2.037419 , 1.5481308 , 1.88460115, 1.40065381,
1.16221919, 1.21493393, 1.1253675 , 1.1257574 , 1.08320837,
2.22150035, 2.63240775, 2.03734889, 1.94949815, 1.41787821,
3.80253946, 3.82716709, 3.02559912, 2.58351668, 3.06317915,
2.51189262, 3.24945661, 3.20053623, 2.31502133, 2.57736614,
4.62272335, 4.70795558, 4.97934873, 5.76980662, 4.97677434,
4.23772352, 4.73154005, 4.31186018, 3.43692443, 2.42336767,
1.89507979, 2.58509517, 1.87114393, 2.09761213, 1.52512991,
1.90355965, 1.14117255, 1.07630129, 1.41866503, 0.87257058,
0.55814673, 2.81302204, 1.98191706, 3.61912165, 2.67814359,
2.84379087, 2.38773647, 1.68431005, 1.48794416, 0.88052901,
1.29770632, 1.73365949, 1.79623326, 1.28510188, 1.06722655,
1.50474152, 1.0926077 , 1.4124716 , 1.27687566, 2.70773448,
3.34342197, 2.70614721, 2.33943225, 2.62586776, 2.02267562,
1.41614464, 1.18733632, 0.76285483, 1.21806093, 2.25246362,
2.45667903, 1.78598427, 2.87602269, 3.75076963, 3.7765598 ,
3.19587481, 3.50146304, 3.53669726, 3.8995017 , 3.06992113,
2.37488444, 2.4049944 , 2.40313759, 2.48375504, 2.58524837,
2.6884848 , 1.98784746, 1.85701545, 2.5028584 , 2.52755525,
2.08079944],
[0.99262075, 1.70123371, 1.42710876, 2.09695269, 1.52332811,
1.38796736, 0.91728086, 1.47248816, 0.97118794, 1.19379776,
1.53992266, 1.11147764, 1.08579688, 0.95727714, 1.09981667,
1.17337987, 1.06915286, 1.16085561, 1.20924113, 1.26272003,
0.82327878, 1.06675543, 1.05250101, 1.19472565, 0.87054942,
1.03161174, 1.47514284, 1.5571645 , 1.16209403, 1.24525565,
1.0997477 , 1.4185143 , 0.89457573, 1.39020645, 2.27128718,
2.53390491, 2.73008011, 1.85626501, 2.63545185, 2.54138971,
2.88738679, 2.96814654, 3.25333939, 2.97501707, 2.62835334,
2.57395103, 2.53805406, 2.72209696, 3.44730891, 3.54841099,
3.27395934, 2.78605488, 2.91952202, 2.38910821, 1.98511209,
4.10186256, 4.4250969 , 4.34459859, 3.78065832, 3.81409911,
3.50859429, 3.33778416, 4.88515458, 4.95815937, 4.14978344,
3.30752754, 2.6854525 , 2.0141847 , 2.17068714, 2.1713234 ,
2.02118407, 1.7895237 , 1.53367833, 3.08507086, 2.32338636,
1.72087692, 2.03807441, 2.73502189, 1.92706757, 1.60367812,
2.88232077, 2.63870686, 2.77199683, 3.46750226, 3.02396473,
2.6863966 , 2.85086389, 2.49723737, 2.24018703, 2.17070273,
1.72242062, 3.14916084, 2.82142156, 1.87669875, 1.93361151,
1.38710944, 1.34383167, 1.24213195, 2.70008131, 2.53972604,
2.30509822],
[0.9991033 , 0.86620851, 1.32749694, 1.11056569, 1.52330344,
1.34240847, 0.97286635, 4.09003129, 3.47848121, 3.13253424,
2.73069016, 2.121885 , 1.68817757, 2.27286062, 2.00810334,
2.18179371, 2.79383138, 2.14213537, 2.94346205, 3.41349862,
4.24659633, 3.64076073, 4.21114023, 3.88988456, 2.99943627,
4.12036555, 3.9215284 , 3.07535748, 2.37102463, 3.87717561,
3.07276223, 2.51940388, 2.69835276, 1.73963556, 1.76738695,
1.25307984, 4.32044268, 5.24881369, 6.15610474, 5.80623441,
5.45150201, 4.98329438, 4.37786238, 4.38188781, 4.0967548 ,
3.34035555, 3.63593765, 2.97775244, 2.97194528, 2.6523777 ,
2.19915773, 2.49332829, 3.37730882, 2.92813687, 2.54077768,
2.27619667, 2.89891564, 2.28779249, 1.84345382, 1.3725854 ,
1.60400574, 1.23894759, 1.36143274, 1.37751795, 1.52542728,
1.1224805 , 1.83308807, 1.86972769, 1.31221991, 1.01753235,
2.00308652, 2.79398193, 3.00997336, 2.86571129, 2.53093516,
2.00761155, 2.53710486, 2.97361304, 2.13461364, 1.40573022,
5.76064129, 4.77818893, 5.49678049, 4.9276452 , 4.84168475,
4.2444369 , 3.88686359, 4.59561311, 4.86865275, 4.45306729,
4.02118307, 3.35321723, 3.44057851, 2.61685051, 2.10466462,
1.85506565, 1.32670471, 2.4786549 , 2.16918857, 2.21181489,
1.87515433]]),
'c': array([[0.88040711, 0.98552687, 0.91932314, 1.05144908, 1.01317215,
1.13608397, 1.2154548 , 1.08726785, 1.17995166, 1.03219313,
0.96747184, 0.98610837, 0.95444347, 0.95458131, 0.93953872,
1.24362044, 1.31041924, 1.21543851, 1.19502874, 1.03636874,
1.51160856, 1.5168326 , 1.38534819, 1.3039992 , 1.39172235,
1.29459402, 1.41618303, 1.40975914, 1.26062355, 1.30319155,
1.603414 , 1.61346078, 1.63922999, 1.73216255, 1.63902196,
1.55068547, 1.61584305, 1.56213896, 1.46333575, 1.28141729,
1.18238609, 1.30420648, 1.17682523, 1.2248495 , 1.07797457,
1.18435616, 0.96003111, 0.93709683, 1.03655949, 0.82027155,
0.55814673, 1.34706239, 1.20256041, 1.48093625, 1.31845181,
1.35358912, 1.27457899, 1.13341938, 1.06295002, 0.82549678,
1.00723617, 1.1448844 , 1.15942175, 1.00418054, 0.93388859,
1.06973683, 0.94286172, 1.03505805, 1.00218631, 1.32472867,
1.43877955, 1.32439198, 1.26530848, 1.30956045, 1.21202957,
1.03594848, 0.97635164, 0.73851383, 0.98721387, 1.24861753,
1.28734375, 1.15704067, 1.35966335, 1.50062708, 1.50609772,
1.40914703, 1.47380543, 1.47594084, 1.52729837, 1.39260766,
1.27211243, 1.27789111, 1.27753475, 1.29089919, 1.30422659,
1.32064541, 1.20393818, 1.17354286, 1.29340771, 1.29665073,
1.22230212],
[0.89909237, 1.13735114, 1.03860645, 1.22474958, 1.07724656,
1.02911763, 0.84962678, 1.05670517, 0.88502032, 0.97863598,
1.08395144, 0.94953291, 0.94045386, 0.87588697, 0.94541035,
0.97141754, 0.93456961, 0.96698977, 0.98409576, 0.99875464,
0.78790825, 0.93372203, 0.92868259, 0.97896403, 0.81894453,
0.92129749, 1.05777777, 1.09091784, 0.96742759, 0.99452085,
0.94538596, 1.03652295, 0.83471938, 1.02966044, 1.25223013,
1.29748452, 1.32946865, 1.17336852, 1.31081897, 1.29846738,
1.36161476, 1.3754826 , 1.41669289, 1.37666238, 1.30988684,
1.30274311, 1.29802936, 1.32777526, 1.46606299, 1.47665076,
1.42053684, 1.34134209, 1.36713293, 1.27484225, 1.20330269,
1.52872815, 1.57815223, 1.56676863, 1.5069671 , 1.5140606 ,
1.47423763, 1.43729891, 1.63135967, 1.63751772, 1.53647294,
1.42935273, 1.32000219, 1.21005693, 1.23592147, 1.23601787,
1.21168304, 1.15786296, 1.08142847, 1.39459702, 1.26222896,
1.14191472, 1.21560706, 1.33051691, 1.1898176 , 1.10971122,
1.36074484, 1.3112464 , 1.33836008, 1.47136629, 1.38506754,
1.32020245, 1.35508945, 1.2926696 , 1.24645176, 1.23592383,
1.14227335, 1.40301287, 1.34884411, 1.17811575, 1.1913379 ,
1.02890965, 1.01841807, 0.99376359, 1.32310527, 1.29824892,
1.25871912],
[0.90334859, 0.81609444, 1.01445813, 0.9492105 , 1.0772366 ,
1.01807305, 0.88612231, 1.52681602, 1.47241259, 1.40082958,
1.32959806, 1.22852721, 1.1343179 , 1.2525321 , 1.20864409,
1.23760428, 1.34299165, 1.23159544, 1.37124384, 1.45718352,
1.55211947, 1.48224771, 1.54638919, 1.52739552, 1.38085558,
1.53171853, 1.52707588, 1.39332152, 1.27137166, 1.52744047,
1.39298073, 1.29558035, 1.32273861, 1.14627278, 1.15272008,
0.99641762, 1.56335264, 1.65803983, 1.79276474, 1.73253051,
1.70703598, 1.63954883, 1.5714726 , 1.57204185, 1.52790265,
1.43797423, 1.4819554 , 1.37713209, 1.3761349 , 1.31304156,
1.24023519, 1.29215628, 1.44767912, 1.36861225, 1.29838701,
1.25317236, 1.36359446, 1.25539782, 1.17039218, 1.02538867,
1.10984359, 0.99299162, 1.02268499, 1.02658444, 1.07809472,
0.95342281, 1.16798397, 1.17649621, 1.01075461, 0.91544847,
1.20747856, 1.34302358, 1.38266498, 1.3578927 , 1.29709456,
1.20852983, 1.29790472, 1.37642129, 1.23045579, 1.03342378,
1.73206997, 1.62055506, 1.72121409, 1.63505193, 1.62696877,
1.55177047, 1.52742603, 1.60012791, 1.62969281, 1.58210764,
1.52606927, 1.44135204, 1.46429541, 1.30837637, 1.22591806,
1.17308988, 1.01426608, 1.29022947, 1.23569441, 1.24215295,
1.17775694]]),
'y': array([[ nan, 1.12534684, 0.78687677, 1.34744014, 0.8937605 ,
1.37130824, 1.44973883, 0.68506838, 1.40069507, 0.66077184,
0.77533547, 1.01044922, 0.88510065, 0.94628717, 0.90347348,
2.07064722, 1.60563385, 0.64926096, 1.08649225, 0.62568533,
3.40195452, 1.42168963, 0.5997479 , 0.8612532 , 1.7196858 ,
0.75686298, 1.97129308, 1.27559897, 0.43470538, 1.47024847,
3.28484003, 1.53768077, 1.73012918, 2.26268195, 0.73724807,
0.73308352, 1.91015011, 1.04037832, 0.54971715, 0.35109955,
0.6960319 , 1.83676679, 0.5262108 , 1.36857751, 0.60872915,
1.43404653, 0.38600889, 0.88610278, 1.27250034, 0.47135976,
0.50323275, 2.81302204, 0.44265944, 2.80079717, 0.43304892,
1.41611451, 0.82302463, 0.51549468, 0.90950896, 0.43428517,
1.23992249, 1.42866582, 1.17801942, 0.61644978, 0.77225915,
1.36473665, 0.63585278, 1.25523832, 0.88059143, 2.41931067,
1.89126587, 0.70627266, 0.88858926, 1.49803779, 0.64055295,
0.56496629, 0.78813035, 0.54132091, 1.19250288, 2.01007421,
1.40264063, 0.55818223, 2.2156319 , 2.15859232, 1.41391012,
0.81188963, 1.62539887, 1.40765677, 1.73570745, 0.57910763,
0.6137053 , 1.24708379, 1.21967914, 1.30187207, 1.33274972,
1.34341194, 0.5516161 , 1.0339107 , 1.78521218, 1.25763202,
0.78834969],
[ nan, 1.60302891, 0.83503207, 1.68902526, 0.60751485,
0.91958174, 0.54048864, 1.40145139, 0.5346158 , 1.10332176,
1.3140028 , 0.63270785, 0.91575492, 0.80466696, 1.014357 ,
1.01125323, 0.85709241, 1.01954319, 1.005682 , 1.02631739,
0.54611512, 1.02961637, 0.91281594, 1.06471631, 0.64399971,
0.9774266 , 1.35931288, 1.11893118, 0.67253503, 1.0408559 ,
0.83647616, 1.25643448, 0.49348481, 1.32735728, 1.89271387,
1.46389501, 1.4318387 , 0.38562299, 1.91841053, 1.15052519,
1.58231834, 1.36608591, 1.58104224, 1.04653824, 0.95008092,
1.18956121, 1.20328573, 1.42007103, 1.98327112, 1.46810277,
1.09861109, 0.83996125, 1.40257359, 0.75909967, 0.81513283,
3.28096268, 1.72330577, 1.35530668, 0.86393687, 1.42672334,
1.09355386, 1.20170966, 2.88964507, 1.54167471, 0.66310971,
0.56355152, 0.71336895, 0.58046187, 1.32635297, 1.18981945,
1.03911327, 0.93954762, 0.87043455, 2.61020851, 0.54838883,
0.60666165, 1.4301641 , 1.87143117, 0.45233734, 0.82956564,
2.36365552, 1.04105213, 1.37816334, 1.96218367, 0.92802197,
0.96555454, 1.41636004, 0.92667422, 0.97539087, 1.12728069,
0.74090278, 2.54000621, 0.9879662 , 0.33049243, 1.20009936,
0.60772215, 0.96772189, 0.90044766, 2.43929453, 1.0939012 ,
1.00154724],
[ nan, 0.76566607, 1.27487716, 0.78187494, 1.35388049,
0.87403829, 0.63231416, 3.99895005, 0.78710518, 1.02616219,
0.91240027, 0.65073829, 0.75015188, 1.69130796, 0.9367584 ,
1.3423615 , 1.80243248, 0.61875365, 1.98739513, 1.7626695 ,
2.19246548, 0.81156001, 1.94470156, 1.09189597, 0.51882278,
2.42085582, 1.20344903, 0.56118234, 0.60488687, 2.72254 ,
0.60554033, 0.75563331, 1.41333805, 0.29524071, 1.14435604,
0.60767962, 4.05094736, 2.35386915, 2.38579219, 1.2247274 ,
1.17411291, 1.05160505, 0.86692956, 1.43517855, 1.14641655,
0.64306079, 1.63843726, 0.71607108, 1.29129392, 0.9767768 ,
0.79285478, 1.48645963, 2.11607821, 0.90202569, 0.90327683,
0.97168647, 1.82474011, 0.67570525, 0.75943942, 0.66587068,
1.23944917, 0.72007733, 1.10317897, 1.02183282, 1.15694709,
0.65278131, 1.6555775 , 1.17136839, 0.58432685, 0.70099379,
1.89589844, 1.95859358, 1.4864671 , 1.15703749, 0.94772565,
0.71207892, 1.69806905, 1.67245289, 0.4575623 , 0.45636448,
5.36971952, 0.54818905, 2.18126492, 0.96330048, 1.38446182,
0.86898513, 1.05956383, 2.11820368, 1.72339329, 1.05215936,
1.00667543, 0.73334774, 1.43312006, 0.54175325, 0.73076677,
0.93238177, 0.61063014, 2.15059433, 0.92134188, 1.23164602,
0.85700929]])}