gsquantify2018.py

import numpy as np


inp = input()
u0,kappa,ita,nu,S,n,K = map(float,inp)
num_sims = S 

t_init = 0.0
t_end  = 10.0
N      = n*10 
dt     = float(t_end - t_init) / N 
y_init = u0

def mu(y, t): 
    return kappa * (ita - y)

def sigma(y, t):
    return nu * (y**0.5)

def sigm_diff(y,t):
    return nu/(2*(y**0.5))
    
def dW(delta_t): 
    """Sample a random number at each call."""
    return np.random.normal(loc = 0.0, scale = np.sqrt(delta_t))

ts    = np.arange(t_init, t_end, dt)
ys    = np.zeros(int(N+1))

ys[0] = y_init

print("P1")
for _ in range(int(num_sims)):
    a = ""
    for i in range(1, ts.size+1):
        t = (i-1) * dt
        y = ys[i-1]
        ys[i] = y + mu(y, t) * dt + sigma(y, t) * dW(dt) + 0.5*sigma(y,t)*sigm_diff(y,t)* (((dW(dt))**2) - dt)
        if i%10 == 1:
            a += str(ys[i]) + " "
    print(a)
print("end")

ts    = np.arange(t_init, t_end, dt)
ys    = np.zeros(int(N+1))

ys[0] = y_init

d = (4*ita*kappa)/(nu**2)

def c(t,s):
    return ((nu**2)*(1 - np.exp(-1*kappa*(t-s))))/(4*kappa)

def lambd(y,t,s):
    return (1/c(t,s))*(np.exp(-1*kappa*(t-s)))*y

def fu(lm,y,z,u):
    if lm + 2*np.log(u) <= 0:
        return 0
    else:
        return (z + (lm + 2*np.log(u) )**0.5)**2 + y**2

print("P2")
a = ""
for i in range(1, int(n)+1):
    y = ys[i-1]
    # if d > 1:
    #     ys[i] = y * c(i*dt,0)*((np.random.normal(lambd(ys,i*dt,0))) + np.random.chisquare(lambd(ys,i*dt,0)-1))
    # else:
    #     ys[i] = y * c(i*dt,0)*(np.random.chisquare(lambd(ys,i*dt,0)) + fu(lambd(ys,i*dt,0),np.random.normal(),np.random.normal(),np.random.uniform()))
    ys[i] = ys[0] * c(i*dt,0) * np.random.noncentral_chisquare(d,lambd(ys[0],i*dt,0))
    a += str(ys[i]) + " "
print(a)
print("end")