# Source code for quantecon.robustlq

"""
Solves robust LQ control problems.

"""
from textwrap import dedent
import numpy as np
from .lqcontrol import LQ
[docs] def robust_rule(self, method='doubling'): """ This method solves the robust control problem by tricking it into a stacked LQ problem, as described in chapter 2 of Hansen- Sargent's text "Robustness." The optimal control with observed state is .. math:: u_t = - F x_t And the value function is :math:-x'Px Parameters ---------- method : str, optional(default='doubling') Solution method used in solving the associated Riccati equation, str in {'doubling', 'qz'}. Returns ------- F : array_like(float, ndim=2) The optimal control matrix from above P : array_like(float, ndim=2) The positive semi-definite matrix defining the value function K : array_like(float, ndim=2) the worst-case shock matrix K, where :math:w_{t+1} = K x_t is the worst case shock """ # == Simplify names == # A, B, C, Q, R = self.A, self.B, self.C, self.Q, self.R beta, theta = self.beta, self.theta k, j = self.k, self.j # == Set up LQ version == # I = np.identity(j) Z = np.zeros((k, j)) if self.pure_forecasting: lq = LQ(-beta*I*theta, R, A, C, beta=beta) # == Solve and convert back to robust problem == # P, f, d = lq.stationary_values(method=method) F = np.zeros((self.k, self.n)) K = -f[:k, :] else: Ba = np.hstack([B, C]) Qa = np.vstack([np.hstack([Q, Z]), np.hstack([Z.T, -beta*I*theta])]) lq = LQ(Qa, R, A, Ba, beta=beta) # == Solve and convert back to robust problem == # P, f, d = lq.stationary_values(method=method) F = f[:k, :] K = -f[k:f.shape[0], :] return F, K, P
[docs] def robust_rule_simple(self, P_init=None, max_iter=80, tol=1e-8): """ A simple algorithm for computing the robust policy F and the corresponding value function P, based around straightforward iteration with the robust Bellman operator. This function is easier to understand but one or two orders of magnitude slower than self.robust_rule(). For more information see the docstring of that method. Parameters ---------- P_init : array_like(float, ndim=2), optional(default=None) The initial guess for the value function matrix. It will be a matrix of zeros if no guess is given max_iter : scalar(int), optional(default=80) The maximum number of iterations that are allowed tol : scalar(float), optional(default=1e-8) The tolerance for convergence Returns ------- F : array_like(float, ndim=2) The optimal control matrix from above P : array_like(float, ndim=2) The positive semi-definite matrix defining the value function K : array_like(float, ndim=2) the worst-case shock matrix K, where :math:w_{t+1} = K x_t is the worst case shock """ # == Simplify names == # A, B, C, Q, R = self.A, self.B, self.C, self.Q, self.R beta, theta = self.beta, self.theta # == Set up loop == # P = np.zeros((self.n, self.n)) if P_init is None else P_init iterate, e = 0, tol + 1 while iterate < max_iter and e > tol: F, new_P = self.b_operator(self.d_operator(P)) e = np.sqrt(np.sum((new_P - P)**2)) iterate += 1 P = new_P I = np.identity(self.j) S1 = P.dot(C) S2 = C.T.dot(S1) K = inv(theta * I - S2).dot(S1.T).dot(A - B.dot(F)) return F, K, P
[docs] def compute_deterministic_entropy(self, F, K, x0): r""" Given K and F, compute the value of deterministic entropy, which is .. math:: \sum_t \beta^t x_t' K'K x_t with .. math:: x_{t+1} = (A - BF + CK) x_t Parameters ---------- F : array_like(float, ndim=2) The policy function, a k x n array K : array_like(float, ndim=2) The worst case matrix, a j x n array x0 : array_like(float, ndim=1) The initial condition for state Returns ------- e : scalar(int) The deterministic entropy """ H0 = np.dot(K.T, K) C0 = np.zeros((self.n, 1)) A0 = self.A - np.dot(self.B, F) + np.dot(self.C, K) e = var_quadratic_sum(A0, C0, H0, self.beta, x0) return e
[docs] def evaluate_F(self, F): """ Given a fixed policy F, with the interpretation :math:u = -F x, this function computes the matrix :math:P_F and constant :math:d_F associated with discounted cost :math:J_F(x) = x' P_F x + d_F` Parameters ---------- F : array_like(float, ndim=2) The policy function, a k x n array Returns ------- P_F : array_like(float, ndim=2) Matrix for discounted cost d_F : scalar(float) Constant for discounted cost K_F : array_like(float, ndim=2) Worst case policy O_F : array_like(float, ndim=2) Matrix for discounted entropy o_F : scalar(float) Constant for discounted entropy """ # == Simplify names == # Q, R, A, B, C = self.Q, self.R, self.A, self.B, self.C beta, theta = self.beta, self.theta # == Solve for policies and costs using agent 2's problem == # K_F, P_F = self.F_to_K(F) I = np.identity(self.j) H = inv(I - C.T.dot(P_F.dot(C)) / theta) d_F = np.log(det(H)) # == Compute O_F and o_F == # AO = np.sqrt(beta) * (A - np.dot(B, F) + np.dot(C, K_F)) O_F = solve_discrete_lyapunov(AO.T, beta * np.dot(K_F.T, K_F)) ho = (np.trace(H - 1) - d_F) / 2.0 tr = np.trace(np.dot(O_F, C.dot(H.dot(C.T)))) o_F = (ho + beta * tr) / (1 - beta) return K_F, P_F, d_F, O_F, o_F