function output = discretion(Omega,A,B,SigmaE,C,Ci,W,beta,n1,n2,X10,Theta1,Theta2);
%-----------------------------------------------------------------------------------------
% PURPOSE: performs discretionary optimization under certainty for Svensson   
%          and Woodford (2002), "Indicator Variables for Optimal Policy"           
%-----------------------------------------------------------------------------------------
% USAGE: output  = discretion(Omega,A,B,SigmaE,C,Ci,W,beta,n1,n2,X10,Theta1,Theta2)
% where:  Omega  = n2xn2 (possibly singular) matrix structural coefficients
%         A      = nxn matrix of structural coefficients for [X(t)' x(t)']
%         B      = nxk matrix of structural coefficients for i(t)
%         SigmaE = n1xn1 covariance matrix of structural shocks e(t)
%         C      = nyxn matrix of target variable coefficients for [X(t)' x(t)']
%         Ci     = nyxk matrix of target variable coefficients for i(t)
%         W      = nyxny matrix of policy weights
%         beta   = discount rate in the central bank's loss function
%         n1     = number of predetermined state variables
%         n2     = number of forward looking state variables
%         X10    = initial value of the state vector X(t) for t = 0
%         Theta1 = ntxn matrix of extra-model coefficients for [X(t)' x(t)']
%         Theta2 = ntxk matrix of extra-model coefficients for i(t)
%-----------------------------------------------------------------------------------------  
%  model setup:
%  [X(t+1)' Omega*x(t+1|t)'] = A*[X(t)' x(t)'] + B*i(t) + [e(t+1)' 0']
%  Y(t) = C*[X(t)' x(t)'] + Ci*i(t)
%  L(t) = Y(t)'*W*Y(t)
%  T(t) = Theta1*[X(t)' x(t)'] + Theta2*i(t)
%
%  model solution:
%  i(t)   = F*X(t)
%  x(t)   = G*X(t)
%  X(t+1) = M*X(t) + e(t+1)
%-----------------------------------------------------------------------------------------
% RETURNS: a structure
% output.F       = F 
% output.G       = G 
% output.M       = M
% output.V       = V
% output.Loss    = Loss 
% output.weights = central bank policy weights
% output.CovT    = SigmaT (Covariance matrix of extra-model variables T(t))
% output.stdT    = stdT (Standard deviations of extra-model varialbes T(t))
%-----------------------------------------------------------------------------------------
% written by: Gregory Givens
%             1 May 2003
%-----------------------------------------------------------------------------------------

n = n1 + n2; N = size(A,1); nn = size(SigmaE,1);
k = size(B,2); h = size(Omega,1); 

%check matrix dimensions%
%----------------------------------------------------------------
if N ~= n;
   error('wrong number of predetermined and forward variables');
end;
if nn ~= n1;
   error('covariance matrix dimensions do not agree');
end;
if h ~= n2;
   error('dimensions of Omega do not agree');
end;
%----------------------------------------------------------------

%construct structural partitioned matrices%
%----------------------------------------------------------------
A_11 = A(1:n1,1:n1); A_12 = A(1:n1,n1+1:n);
A_21 = A(n1+1:n,1:n1); A_22 = A(n1+1:n,n1+1:n);
B_1 = B(1:n1,:); B_2 = B(n1+1:n,:);
C_1 = C(:,1:n1); C_2 = C(:,n1+1:n);
%----------------------------------------------------------------

%construct partitioned matrix of policy weights%
%----------------------------------------------------------------
Q = C'*W*C; U = C'*W*Ci; R = Ci'*W*Ci;
Q_11 = Q(1:n1,1:n1); Q_12 = Q(1:n1,n1+1:n);
Q_21 = Q(n1+1:n,1:n1); Q_22 = Q(n1+1:n,n1+1:n);
U_1 = U(1:n1,:); U_2 = U(n1+1:n,:);
%----------------------------------------------------------------

%make initial guess of (F(t+1),G(t+1),V(t+1))%
%----------------------------------------------------------------
Ft = zeros(k,n1);
Gt = zeros(n2,n1);
Vt = eye(n1)*.1;
%----------------------------------------------------------------

%solve for the fixed point (F,G,V) via t -> -inf%
%----------------------------------------------------------------
dd = 1; tol = 10^(-8); count = 0;
while (dd>tol)&(count<10000);
   Fto = Ft; Gto = Gt; Vto = Vt;
   Atil = inv(A_22-Omega*Gt*A_12)*(Omega*Gt*A_11-A_21);
   Btil = inv(A_22-Omega*Gt*A_12)*(Omega*Gt*B_1-B_2);
   Astar = A_11 + A_12*Atil; Bstar = A_12*Btil + B_1;
   Qstar = Q_11 + Q_12*Atil + Atil'*Q_21 + Atil'*Q_22*Atil;
   Ustar = Q_12*Btil + Atil'*Q_22*Btil + U_1 + Atil'*U_2;
   Rstar = R + Btil'*Q_22*Btil + Btil'*U_2 + U_2'*Btil;
   Ft = -inv(Rstar+beta*Bstar'*Vt*Bstar)*(Ustar'+beta*Bstar'*Vt*Astar);
   Gt = Atil + Btil*Ft;
   Vt = Qstar+Ustar*Ft+Ft'*Ustar'+Ft'*Rstar*Ft+beta*(Astar+Bstar*Ft)'*Vt*(Astar+Bstar*Ft);
   count = count+1;
   X = max(abs(Fto-Ft)); Y = max(abs(Gto-Gt)); Z = max(abs(Vto-Vt));
   XX = norm(X); YY = norm(Y); ZZ = norm(Z); 
   dd = XX + YY + ZZ;
end;
if count >= 10000;
   disp('max number of optimization iterations reached');
end;
M = Astar + Bstar*Ft;
%----------------------------------------------------------------

%calculate loss%
%----------------------------------------------------------------
delta = (1/(1-beta))*beta*trace(Vt*SigmaE);
Loss = X10'*Vt*X10 + delta;
%----------------------------------------------------------------

%calculate 2nd moments for extra-model variables%
%----------------------------------------------------------------
Theta1_1 = Theta1(:,1:n1); Theta1_2 = Theta1(:,n1+1:n);
Thetah = Theta1_1 + Theta1_2*Gt + Theta2*Ft;
vec_SigmaX = inv(eye(n1^2)-kron(M,M))*vec(SigmaE);
SigmaX = reshape(vec_SigmaX,n1,n1);
SigmaT = Thetah*SigmaX*Thetah';
stdT = sqrt(diag(SigmaT)');
%----------------------------------------------------------------

%input results into 'output' structure%
%----------------------------------------------------------------
output.meth = 'optimal discretion: complete information';
output.F = Ft;
output.G = Gt;
output.M = M;
output.V = Vt;
output.Loss = Loss; output.weights = diag(W)';
output.CovT = SigmaT; output.stdT = stdT;
%----------------------------------------------------------------

Ag = [eye(n1) zeros(n1,n2);zeros(n2,n1) Omega];
Bg = A + B*[Ft zeros(k,n2)];
egv = abs(eig(Bg,Ag));
output.numb = sum((egv>=1));