Consider the risk-neutral interest rate process drt = a(b rt) dt + rt dZt (1) where a, b, are positive constants and Zt is a standard Brownian motion....

Consider the risk-neutral interest rate process drt = a(b − rt) dt + σ √ rt dZt (1) where a, b, σ are positive constants and Zt is a standard Brownian motion. A European option written on the interest rate rt, with payoff P(r(T)) at the expiry T can be valued by determining 

option value = E Q [exp (-∫ r(t) dt) P(r(T))] . (2) where E Q(·) denotes expectation assuming rt following (1). 

Provide a precise pseudo (NOT Matlab) code to compute the value of the option starting at r(0) = r0, t = 0 using a Monte Carlo method with M simulation paths and N time steps. You may assume the existence of a function which returns random sample φ for a standard normal and a function to calculate the payoff.