Time Dependent Condition not working in Custom PDE

[AnttiZhang0]

Hi
I have a time dependent Boundary Condition, which is Time and X dependent, so it code below py function

class BlackScholesForwardPDE(PDEBase):

  def BC_EuropeanCall_up(self, adjacent_value, dx, x, t):
          Spot = np.exp(x) * self.terms_dict['Spot0']
          BC_state = (Spot - self.terms_dict['Strike']) * np.exp(-1 * 0.05 * (t))
          return BC_state
  def set_BoundaryCondition(self,Option='EuropeanCall'):
        if Option=='EuropeanCall':
            self.bc = [{"value_expression": self.BC_EuropeanCall_down}, {"value_expression": self.BC_EuropeanCall_up}]

Also for a Custom PDE, i also code the evolution_rate

    def evolution_rate(self, state,t):
        gradient_x = state.gradient(bc=self.bc,)[0]
        laplace_x = state.laplace(bc=self.bc)
        x = state.grid.cell_coords[:, 0]
        volatility = self.get_volatility(x,t)  # a local function of x and t
        rate = self.get_r(x,t) # a local function of x and t
        dividend = self.get_q(x,t) # a local function of x and t
        dvdtao = 0.5 * volatility**2 * laplace_x + (rate -dividend - 0.5 *  volatility**2) * gradient_x - rate * state
        return dvdtao

What i am expecting is that: when the py-pde is solving the PDE, it will call the evolution_rate at different time t. At this given time t, when conduct the state.gradient operation, the Boundary Function shall be evaluate at that t. But when the evolution_rate function is called with any t<>0, the BC_EuropeanCall_up always get called as t=0, not matter the t used to evaluate the evolution_rate
it seems below code did not capture t in *args

                def virtual_from_value(adjacent_value, *args):
                    return 2 * value_func(adjacent_value, *args) - adjacent_value

Full code

class FXO_Payoffs:

    def __init__(self,Terms_dict):
        self.Terms_dict =Terms_dict

    def create_InitialState(self,x_array,Option):
        InitialState = np.zeros(x_array.shape)
        for i in range(0,x_array.shape[0]):
            x = x_array[i]
            Spot = np.exp(x) * Terms_dict['Spot0']
            Payoff = self.get_Payoff(Option,Spot)
            InitialState[i] = Payoff
        return InitialState
        
    def get_Payoff(self,Option,Spot):
        if Option=='EuropeanCall':
            Payoff = self.EuropeanCall_Payoff(Spot=Spot)
        return Payoff
    
    def EuropeanCall_Payoff(self,Spot):
        Payoff = max(0,Spot-self.Terms_dict['Strike']) * self.Terms_dict['Notional']
        return Payoff

    def EuropeanPut_Payoff(self,Spot):
        Payoff = max(0,self.Terms_dict['Strike']-Spot)* self.Terms_dict['Notional']
        return Payoff

class BlackScholesForwardPDE(PDEBase):

    def __init__(self,Vol_array,Underlying_array,Time_array,Rates_array,Terms_dict):
        self.Vol_array = Vol_array  # as In(S) X t(inverse time)
        self.Underlying_array = Underlying_array
        self.Time_array = Time_array
        self.Rates_array = Rates_array
        self.Terms_dict = Terms_dict

    def set_BoundaryCondition(self,Option='EuropeanCall'):
        if Option=='EuropeanCall':
            self.bc = [{"value_expression": self.BC_EuropeanCall_down}, {"value_expression": self.BC_EuropeanCall_up}]

    def BC_EuropeanCall_up(self, adjacent_value, dx, x, t):
        Spot = np.exp(x) * self.Terms_dict['Spot0']
        Rate = self.get_r(x_array=x,t=t)
        BC_state = (Spot - self.Terms_dict['Strike']) * np.exp(-1 * Rate * (t)) * self.Terms_dict['Notional']
        return BC_state

    def BC_EuropeanCall_down(self, adjacent_value, dx, x, t):
        BC_state = 0* self.Terms_dict['Notional']
        return BC_state

    def get_r(self,x_array,t):
        Rate_array = np.ones((x_array.shape))
        for j in range(0,len(self.Time_array)-1):
            if self.Time_array[j] <= t and self.Time_array[j+1] >= t:
                Rate_short = self.Rates_array[j]
                Rate_long = self.Rates_array[j]
                Time_short = self.Time_array[j]
                Time_long = self.Time_array[j + 1]
                break
        Rate = Rate_short + (Rate_long - Rate_short) * (t-Time_short) * (Time_long-Time_short)
        return Rate_array * Rate

    def get_q(self,x_array,t):
        Dividend_array = np.ones((x_array.shape))
        for j in range(0,len(self.Time_array)-1):
            if self.Time_array[j] <= t and self.Time_array[j+1] >= t:
                Dividend_short = self.Rates_array[j]
                Dividend_long = self.Rates_array[j]
                Time_short = self.Time_array[j]
                Time_long = self.Time_array[j + 1]
                break
        Dividend = Dividend_short + (Dividend_long - Dividend_short) * (t-Time_short)* (Time_long-Time_short)
        return Dividend_array * Dividend

    def get_volatility(self,x_array,t):
        vol_array = np.zeros((x_array.shape))
        for k in range(0,x_array.shape[0]):
            x = x_array[k]
            for i in range(0,len(self.Underlying_array)-1):
                if self.Underlying_array[i] <= x and self.Underlying_array[i+1] >= x:
                    Underlying_short = self.Underlying_array[i]
                    Underlying_long = self.Underlying_array[i+1]
                    break
            for j in range(0,len(self.Time_array)-1):
                if self.Time_array[j] <= t and self.Time_array[j+1] >= t:
                    Time_short = self.Time_array[j]
                    Time_long = self.Time_array[j+1]
                    break
            Vol_x0t0 = self.Vol_array[i,j]
            Vol_x1t0 = self.Vol_array[i+1, j]
            Vol_x0t1 = self.Vol_array[i, j+1]
            Vol_x1t1 = self.Vol_array[i, j]
            Vol_t0 = Vol_x0t0 + (Vol_x1t0 - Vol_x0t0)* (x-Underlying_short)/(Underlying_long-Underlying_short)
            Vol_t1 = Vol_x0t1 + (Vol_x1t1 - Vol_x0t1) * (x - Underlying_short) / (Underlying_long - Underlying_short)
            Vol_x = Vol_t0 + (Vol_t1 - Vol_t0) * (t-Time_short) * (Time_long-Time_short)
            vol_array[k] = Vol_x
        return vol_array

    def evolution_rate(self, state,t):
        # translate from backward to forward, by using negative t, t = backward time from T, x is the S, state is V
        gradient_x = state.gradient(bc=self.bc,t=t)[0]
        laplace_x = state.laplace(bc=self.bc)
        x = state.grid.cell_coords[:, 0]
        volatility = self.get_volatility(x,t)
        rate = self.get_r(x,t)
        dividend = self.get_q(x,t)
        dvdtao = 0.5 * volatility**2 * laplace_x + (rate -dividend - 0.5 *  volatility**2) * gradient_x - rate * state
        return dvdtao

# define grid and an initial state
support_points = 100
x_min = np.log(0.5)
x_max = np.log(2.0)
PDE_grid = CartesianGrid(bounds=[(x_min,x_max)], shape=[support_points])
field = ScalarField(PDE_grid)
Vol_array = np.ones((support_points,support_points)) * 0.08
Underlying_array = np.squeeze(PDE_grid.cell_coords)
Time_array = np.array(np.linspace(start=0,stop=2, num=support_points))
Rates_array = np.array(np.linspace(start=0.05,stop=0.05, num=support_points))
Terms_dict = {
    'Spot0': 1.1,
    'Strike': 1.1,
    'T':2.0,
    'Notional': 1000000,
}
BS_PDE = BlackScholesForwardPDE(Vol_array=Vol_array,Underlying_array=Underlying_array,Time_array=Time_array,Rates_array=Rates_array,Terms_dict=Terms_dict)
BS_PDE.set_BoundaryCondition()
#y = BS_PDE.evolution_rate(state=field,t=1)
FXO_Payoffs0 = FXO_Payoffs(Terms_dict=Terms_dict)
FinalPayoff = FXO_Payoffs0.create_InitialState(x_array=Underlying_array,Option='EuropeanCall')
storage = MemoryStorage()
InitialState = ScalarField(grid=PDE_grid,data=FinalPayoff)
BS_PDE.solve(state=InitialState, t_range=(0, 1.0), tracker=storage.tracker(0.1))
plot_kymograph(storage)

[david-zwicker]

I must admit I don't understand your question. What output do you expect and what do you get instead?

[AnttiZhang0]

Hi David
Apologize for my expression, I have edit the question to make it easier to understand. Thank you

[david-zwicker]

If I understand correctly, you want to change the boundary conditions as a function of time, but I don't see where in evolution_rate you are calling any function that sets a time-dependent boundary condition. I don't have the time to look at your code in detail, but if you have a minimal example, I might be able to see the problem better.