function [phi_star,pi_star,pi_star1,pi_star2,phi_og,phi_within,phi_across,phi_adjust,phi_dc,phi_dc_eff] = optimalpolicy(mu, theta, sigma_c, sigma, varepsilon, Omega, alpha, beta, gamma, varphi)

    N = length(Omega);
    
    beta                  = [beta;0];

    % Cost-based IO variables

    Omega_tilde           = [Omega, alpha; 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                    = [mu;1];
    Omega                 = diag(mu)^(-1) * Omega_tilde;
    Psi                   = (eye(N+1)-Omega)^(-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(j) * weightedcov(Psi_tilde(:,k), Psi(:,N+1)/Lambda, Omega_tilde(j,:));
        end
    end

    J        = M .* (1-theta.^(-1));
    phi_dc   = Domar_tilde'* (diag(theta)^(-1)-eye(N)) + varphi * J';
    phi_dc_eff = Domar_tilde'* (diag(theta)^(-1)-eye(N));
    xi       = 1 - Lambda_hat - J' * Psi_hat(1:N,N+1);
    H1       = diag(varepsilon) * diag(Domar.* Psi(1:N,N+1)/Lambda - M) * (diag(theta)^(-1)-eye(N));
    B        = zeros(N,N);

    for k = 1:N
        for l = 1:N
            B(k,l) = sigma_c * weightedcov(Psi(:,k), Psi(:,l), beta);
            for j = 1:N
                B(k,l) = B(k,l) + sigma(j) * Domar(j) * weightedcov(Psi(:,k), Psi(:,l), Omega(j,:));
            end
        end
    end

    H2 = (diag(theta)^(-1)-eye(N)) * B * (diag(theta)^(-1)-eye(N));
    phi_og       = Domar_tilde'* (diag(theta)^(-1)-eye(N)) - J';
    phi_within   = (gamma + varphi)/xi.* (Psi_hat(1:N,N+1))' * H1;
    phi_across   = (gamma + varphi)/xi.* (Psi_hat(1:N,N+1))' * H2;
    phi_adjust   = (1+varphi)*(gamma-1)*(1-Lambda_hat-xi)/xi * J';
    phi_star     = phi_og + phi_within + phi_across + phi_adjust;
    pi_star1     = 1 - Lambda;
    pi_star2     = (1 + varphi + (gamma -1) * Lambda) * ((1-Lambda_hat)/xi -1);
    pi_star      = pi_star1 + pi_star2;



    
end