function [v] = BS_Theta(PC, S, X, T, v, r, b)

[d1, d2] = BS_D1D2(S, X, T, v);

v = -(S * exp(-((PC * d1) ^ 2) / 2) / sqrt(2 * pi) * v * exp(- r * T)) / (2 * sqrt(T)) + PC * r * S * BS_CND(PC * d1) * exp(- r * T) - PC * r * X * exp(-r * T) * BS_CND(PC * d2);
