% This function is used to solve for the output gap given realized
% productivity shocks, the policy weights and the inflation bias
function y_tilde = solveoutputgap(phi_star, pi_star, dlogz)

global N mu theta sigma_c sigma omega_tilde alpha_tilde beta_tilde gamma varphi

beta                  = [beta_tilde;0];

% Cost-based IO variables

Omega_tilde           = [omega_tilde, alpha_tilde; zeros(1,N+1)];
Psi_tilde             = (eye(N+1)-Omega_tilde)^(-1);
Domar_tilde           = (beta' * Psi_tilde)';
Lambda_tilde          = Domar_tilde(N+1);
Domar_tilde           = Domar_tilde(1:N);

% Revenue-based IO variables

mu_trans              = [mu;1];
Psi                   = (eye(N+1)-diag(mu_trans)^(-1) * Omega_tilde)^(-1);
Domar                 = (beta' * Psi)';
Lambda                = Domar(N+1);
Domar                 = Domar(1:N);

% Sticky-price IO variables

theta_trans           = [theta;1];
Psi_hat               = (eye(N+1)-diag(theta_trans) * Omega_tilde)^(-1);
Domar_hat             = (beta' * Psi_hat)';
Lambda_hat            = Domar_hat(N+1);
Domar_hat             = Domar_hat(1:N);


M        = zeros(N,1);

for k = 1:N
    M(k) = sigma_c * weightedcov(Psi_tilde(:,k), Psi(:,N+1)/Lambda, beta);
    for j = 1:N
        M(k) = M(k) + sigma(j) * Domar(j)/mu_trans(j) * weightedcov(Psi_tilde(:,k), Psi(:,N+1)/Lambda, Omega_tilde(j,:));
    end
end

J        = M .* (1-theta.^(-1));
xi       = 1 - Lambda_hat - J' * Psi_hat(1:N,N+1);
K        = (gamma + varphi)/(1-Lambda_hat + varphi * (1-Lambda_hat-xi)) * Psi_hat(1:N,N+1);
V        = 1/(1-Lambda_hat + varphi * (1-Lambda_hat-xi)) * Psi_hat(1:N,N+1) * (Domar_tilde' - Domar_hat' * diag(theta) + varphi*J'*Psi_hat(1:N,1:N)*diag(theta)) - Psi_hat(1:N,1:N)*diag(theta);
y_tilde  = (pi_star - phi_star * V * dlogz)/(phi_star * K);
% dlogw    = 1/(1-Lambda_hat + varphi * (1 - Lambda_hat - xi)) * ((gamma + varphi) * y_tilde + (Domar_tilde' - Domar_hat' * diag(theta) + varphi*J'*Psi_hat(1:N,1:N)*diag(theta)) * dlogz);

end

