diff --git a/models/BJ_22/BJ_22.json b/models/BJ_22/BJ_22.json
new file mode 100644
index 0000000..123d0db
--- /dev/null
+++ b/models/BJ_22/BJ_22.json
@@ -0,0 +1,32 @@
+{
+  "name": "BJ_22",
+  "author": "Brotherhood and Jerbashian",
+  "title": "Firm behavior during an epidemic",
+  "year": 2022,
+  "publication": "Under Review at Journal of Economic Dynamics and Control",
+  "variabledim": 1,
+  "initial_type": "Param",
+  "initial_name": "initial_infected",
+  "initial_value": 0.001,
+  "Code_type": "Matlab",
+  "Dynare_version": "0",
+  "model_type":2, 
+  "features":{
+		"fiscal_policy": 1, 
+		"monetary_policy": 0, 
+		"financial": 0, 
+		"firm_exit": 0, 
+		"multi_country": 0, 
+		"heterogeneity": 0, 
+		"optimal_policy": 1 
+  }, 
+
+  "variables": [ 
+			"Labour",
+			"Output",
+			"Susceptibles",
+			"Infected",
+			"Recovered",
+			"Deaths"
+			]
+}
diff --git a/models/BJ_22/compute_profits.m b/models/BJ_22/compute_profits.m
new file mode 100644
index 0000000..4c25dcb
--- /dev/null
+++ b/models/BJ_22/compute_profits.m
@@ -0,0 +1,555 @@
+function [out, outstruct] = compute_profits(in)
+    
+    % "in" is T by 2 (no inner concavity) or T by 3 (with inner concavity)
+
+    global N T q q_m1 rho_r_s rho_r_a w delta_n delta_h delta_s delta_l ...
+        A beta alpha1 alpha2 gamma rho_d initial_infected varphi Pi_pq Pi_pn ...
+        h_over_n_plus_h n_nd m_nd flag_internalize_q Pi_q flag_planner ...
+        delta_d planner_obj_n
+    
+
+    % pre-allocates some vectors
+    n_vec = zeros(T, 1);
+%     m_vec = zeros(T, 1);
+%     z_vec = zeros(T, 1);
+    
+    %---------------------------------------------------------------------%
+    %             Distribution of workers across health states            %
+    %---------------------------------------------------------------------%    
+
+    % pre-allocates some variables
+    if nargout>1
+        n_vec = zeros(T, 1);
+        l_vec = zeros(T, 1);
+        a_vec = zeros(T, 1);
+        s_vec = zeros(T, 1);
+        r_s_vec = zeros(T, 1);
+        r_a_vec = zeros(T, 1);
+        d_vec = zeros(T, 1);
+        nn_vec = zeros(T, 1);
+        mn_vec = zeros(T, 1);
+        nm_vec = zeros(T, 1);
+        mm_vec = zeros(T, 1);
+        nr_vec = zeros(T, 1);
+        mr_vec = zeros(T, 1);
+        production_vec = zeros(T, 1);
+        profit_vec = zeros(T, 1);
+        h_vec = zeros(T, 1);
+        infected_workplace_vec = zeros(T, 1);
+        p_vec = zeros(T, 1);
+        n_tilde_vec = zeros(T, 1);
+        m_tilde_vec = zeros(T, 1);
+    end
+    
+    % starts profit
+    disc_profits = 0;
+    disc_planner_obj = 0;
+    
+    % loop to compute profits
+    for t = 1:T
+        % computes state variables
+        if t==1
+            i_tm1 = initial_infected * N; % initial conditions
+            i_n_tm1 = i_tm1 * n_nd / (n_nd + m_nd);
+            i_m_tm1 = i_tm1 * m_nd / (n_nd + m_nd);
+            z_tm1 = 0;
+            d_tm1 = 0;
+            a_tm1 = 0;
+            a_n_tm1 = 0;
+            a_m_tm1 = 0;
+            r_s_tm1 = 0;
+            r_a_tm1 = 0;
+            r_a_n_tm1 = 0;
+            r_a_m_tm1 = 0;
+            nn_tm1 = n_nd; % I think I don't need to store all of this...
+            nm_tm1 = 0;
+            mn_tm1 = 0;
+            mm_tm1 = m_nd;%until here
+            q_tm1 = q_m1; % check this (use formula for q)            
+            
+%             n_tm1 = n_nd; % initial conditions
+%             m_tm1 = m_nd;
+%             n_tilde_tm1 = (1 - h_over_n_plus_h) * initial_infected;
+%             m_tilde_tm1 = h_over_n_plus_h * initial_infected;
+%             s_tm1 = 0;
+%             a_tm1 = 0;
+%             r_s_tm1 = 0;
+%             r_a_tm1 = 0;
+%             d_tm1 = 0;
+        else
+            i_tm1 = i_t;
+            i_n_tm1 = i_n_t;
+            i_m_tm1 = i_m_t;
+            z_tm1 = z_t;
+            d_tm1 = d_t;
+            a_tm1 = a_t;
+            a_n_tm1 = a_n_t;
+            a_m_tm1 = a_m_t;
+            r_s_tm1 = r_s_t;
+            r_a_tm1 = r_a_t;
+            r_a_n_tm1 = r_a_n_t;
+            r_a_m_tm1 = r_a_m_t;
+            nn_tm1 = nn_t; % I think I don't need to store all of this...
+            nm_tm1 = nm_t;
+            mn_tm1 = mn_t;
+            mm_tm1 = mm_t; % ... until here
+            q_tm1 = q(t-1);               
+            
+            
+            
+%             n_tm1 = n_t;
+%             m_tm1 = m_t;
+%             n_tilde_tm1 = n_tilde_t;
+%             m_tilde_tm1 = m_tilde_t;
+%             s_tm1 = s_t;
+%             a_tm1 = a_t;
+%             r_s_tm1 = r_s_t;
+%             r_a_tm1 = r_a_t;
+%             d_tm1 = d_t;
+        end
+        
+%         x_nn_vec = in(t, 1);
+%         x_nm_vec = in(t, 2);        
+        
+        x_nn_t = in(t, 1);
+%         x_nn_t = 0.5 + 0.5 * sin(in(t, 1));
+        x_mn_t = 1 - x_nn_t;
+
+        x_nm_t = in(t, 2);
+%         x_nm_t = 0.5 + 0.5 * sin(in(t, 2));
+        x_mm_t = 1 - x_nm_t;
+
+        x_nr_t = 1;
+        x_mr_t = 1 - x_nr_t;
+
+        % choices in levels
+        nn_t = x_nn_t * (nn_tm1 + nm_tm1 - varphi * i_n_tm1);
+        mn_t = x_mn_t * (nn_tm1 + nm_tm1 - varphi * i_n_tm1);
+
+        nm_t = x_nm_t * (mn_tm1 + mm_tm1 - varphi * i_m_tm1);
+        mm_t = x_mm_t * (mn_tm1 + mm_tm1 - varphi * i_m_tm1);
+        
+        if isnan(nn_t)
+            disp('stop here')
+        end
+        
+        if t < T
+            r_s_t = r_s_tm1 + (1-rho_d) * rho_r_s * z_tm1;
+        else
+            r_s_t = r_s_tm1 + z_tm1; % cure
+        end
+
+        nr_t = x_nr_t * r_s_t;
+        mr_t = x_mr_t * r_s_t;     
+        
+%         r_s_t = r_s_t_hat - fr_t + h_t * k_r_s_t;
+
+        nr_t = x_nr_t * r_s_t;
+        mr_t = x_mr_t * r_s_t;
+
+        % infection probability of on-site worker
+        p_tm1 = min(Pi_pq * q_tm1 + Pi_pn * (i_n_tm1 + a_n_tm1), 1);
+
+        % susceptible and exposed in t-1
+        s_n_tm1 = nn_tm1 + nm_tm1 - i_n_tm1 - a_n_tm1 - r_a_n_tm1;
+        e_n_tm1 = s_n_tm1 * p_tm1;
+
+        if s_n_tm1 < 0
+            disp('erro')
+        end
+        
+        s_m_tm1 = mn_tm1 + mm_tm1 - i_m_tm1 - a_m_tm1 - r_a_m_tm1;
+        e_m_tm1 = s_m_tm1 * q_tm1;
+        
+        if s_m_tm1 < 0
+            disp('erro')
+        end
+
+        % uncertain in t before hiring/firing
+        s_hat_n_t = s_n_tm1 * (1-p_tm1);
+        i_hat_n_t = e_n_tm1;
+        a_hat_n_t = (1-rho_r_a) * a_n_tm1 + (1-varphi) * i_n_tm1;
+        r_a_hat_n_t = r_a_n_tm1 + rho_r_a * a_n_tm1;
+
+        s_hat_m_t = s_m_tm1 * (1-q_tm1);
+        i_hat_m_t = e_m_tm1;
+        a_hat_m_t = (1-rho_r_a) * a_m_tm1 + (1-varphi) * i_m_tm1;
+        r_a_hat_m_t = r_a_m_tm1 + rho_r_a * a_m_tm1;
+        
+        % some measures of workers by health states in t
+        d_t = d_tm1 + rho_d * z_tm1;
+        r_a_t = r_a_tm1 + rho_r_a * a_tm1;% - f_t_r_a + h_t * k_r_a_t;
+        a_t = (1-rho_r_a) * a_tm1 + (1-varphi) * i_tm1;% - f_t_a + h_t * k_a_t;
+        
+        if t < T
+            z_t = (1-rho_d) * (1-rho_r_s) * z_tm1 + varphi * i_tm1;
+        else
+            z_t = 0; % cure
+        end
+        
+        if z_t < 0
+            disp('stop here')
+        end
+
+        % incubated in t
+        if nn_t + mn_t < 1e-10
+            prop_nn_t = 0;
+            prop_mn_t = 0;
+        else
+            prop_nn_t = nn_t / (nn_t + mn_t);
+            prop_mn_t = mn_t / (nn_t + mn_t);
+        end
+        
+        if nm_t + mm_t < 1e-10
+            prop_nm_t = 0;
+            prop_mm_t = 0;
+        else
+            prop_nm_t = nm_t / (nm_t + mm_t);
+            prop_mm_t = mm_t / (nm_t + mm_t);
+        end
+        
+%         if nh_t + mh_t < 1e-10
+            prop_nh_t = 0;
+            prop_mh_t = 0;
+%         else
+%             prop_nh_t = nh_t / (nh_t + mh_t);
+%             prop_mh_t = mh_t / (nh_t + mh_t);
+%         end
+        
+        i_nn_t = (e_n_tm1) * prop_nn_t;
+        i_nm_t = (e_m_tm1) * prop_nm_t;
+
+        i_n_t = i_nn_t + i_nm_t;
+        
+        if isnan(i_n_t)
+            disp('stop here')
+        end
+
+        i_mn_t = (e_n_tm1) * prop_mn_t;
+        i_mm_t = (e_m_tm1) * prop_mm_t;
+
+        i_m_t = i_mn_t + i_mm_t;
+
+        i_t = i_n_t + i_m_t;
+        
+        if i_t < 0
+            disp('stop here')
+        end
+
+        % asymptomatics
+%         if nn_t + mn_t < 1e-10
+%             prop_nn_t = 0;
+%             prop_mn_t = 0;
+%         else
+%             prop_nn_t = nn_t / (nn_t + mn_t);
+%             prop_mn_t = mn_t / (nn_t + mn_t);
+%         end
+%         
+%         if nm_t + mm_t < 1e-10
+%             prop_nm_t = 0;
+%             prop_mm_t = 0;
+%         else
+%             prop_nm_t = nm_t / (nm_t + mm_t);
+%             prop_mm_t = mm_t / (nm_t + mm_t);
+%         end
+        
+        a_nn_t = (a_hat_n_t) * prop_nn_t;
+        a_mn_t = (a_hat_n_t) * prop_mn_t;
+        a_nm_t = (a_hat_m_t) * prop_nm_t;
+        a_mm_t = (a_hat_m_t) * prop_mm_t;
+
+        a_n_t = a_nn_t + a_nm_t;
+        a_m_t = a_mn_t + a_mm_t;
+
+        % recovered asymptomatics
+        r_a_nn_t = (r_a_hat_n_t) * prop_nn_t;
+        r_a_mn_t = (r_a_hat_n_t) * prop_mn_t;
+        r_a_nm_t = (r_a_hat_m_t) * prop_nm_t;
+        r_a_mm_t = (r_a_hat_m_t) * prop_mm_t;
+
+        r_a_n_t = r_a_nn_t + r_a_nm_t;
+        r_a_m_t = r_a_mn_t + r_a_mm_t;       
+        
+        n_t = nn_t + nm_t;
+        m_t = mn_t + mm_t;        
+        
+        n_tilde_t = i_nn_t + i_mn_t;
+        m_tilde_t = i_nm_t + i_mm_t;
+        
+        
+        % choices
+        n = nn_t + nm_t + nr_t;
+        m = mn_t + mm_t + mr_t;        
+        
+        if n < 0
+            n = 0;
+        end
+        
+%         % p and q
+%         if t==1
+%             q_tm1 = q_m1;
+%         else
+%             q_tm1 = q(t-1);
+%         end        
+        
+%         % some state variables
+%         r_s_t = r_s_tm1 + (1 - rho_d) * rho_r_s * s_tm1;
+%         r_a_t = r_a_tm1 + rho_r_a * a_tm1;
+%         s_t = (1 - rho_d) * (1 - rho_r_s) * s_tm1 + varphi * (n_tilde_tm1 + m_tilde_tm1);
+%         a_t = (1 - rho_r_a) * a_tm1 + (1 - varphi) * (n_tilde_tm1 + m_tilde_tm1);
+%         d_t = d_tm1 + rho_d * s_tm1;
+%         
+%         % choices
+%         xnn_t = in(t, 1);
+%         xnm_t = in(t, 2);
+% %         xnn_t = sin(in(t, 1)) * 0.5 + 0.5;
+% %         xnm_t = sin(in(t, 2)) * 0.5 + 0.5;
+%         
+%         nn_t =      xnn_t  * (n_tm1 - varphi * n_tilde_tm1);
+%         mn_t = (1 - xnn_t) * (n_tm1 - varphi * n_tilde_tm1);
+%         nm_t =      xnm_t  * (m_tm1 - varphi * m_tilde_tm1);
+%         mm_t = (1 - xnm_t) * (m_tm1 - varphi * m_tilde_tm1);
+%         
+%         if alpha1==1
+%             xnr_t = 1;
+%         else
+%             xnr_t = in(t, 3);
+% %             xnr_t = sin(in(t, 3)) * 0.5 + 0.5;
+%         end
+%         
+%         nr_t =      xnr_t  * r_s_t;
+%         mr_t = (1 - xnr_t) * r_s_t;
+%         
+%         % p and q
+%         if t==1
+%             q_tm1 = q_m1;
+%         else
+%             q_tm1 = q(t-1);
+%         end
+%         
+%         a_n_tm1 = n_tm1 * (a_tm1/(N - r_s_tm1 - s_tm1 - d_tm1));
+% %         a_m_tm1 = m_tm1 * (a_tm1/(N - r_s_tm1 - s_tm1 - d_tm1));
+%         
+%         p_tm1 = Pi_pq * q_tm1 + Pi_pn * (a_n_tm1 + n_tilde_tm1);
+%         p_tm1 = min(p_tm1, 1);
+%         
+%         % other state variables
+%         n_t = nn_t + nm_t;
+%         m_t = mn_t + mm_t;
+%         n_tilde_t = (nn_t * p_tm1 + nm_t * q_tm1) * (1 - (r_a_t + a_t)/(N - r_s_t - s_t - d_t));
+%         m_tilde_t = (mn_t * p_tm1 + mm_t * q_tm1) * (1 - (r_a_t + a_t)/(N - r_s_t - s_t - d_t));
+        
+        % prop. of infected in the workplace
+        if n_tilde_t + m_tilde_t > 1e-6
+            infected_workplace = (nn_t + mn_t) * (p_tm1 - q_tm1 * Pi_pq) * ...
+                (1 - (r_a_t + a_t)/(N - r_s_t - z_t - d_t));
+%         if p_tm1 ~= q_tm1
+%             zz = 1;
+%         end
+            infected_workplace = infected_workplace/(n_tilde_t + m_tilde_t);
+        else
+            infected_workplace = NaN;
+        end
+%         
+%         % choices
+%         n = nn_t + nm_t + nr_t;
+%         m = mn_t + mm_t + mr_t;
+        
+        % splitting m among h and l
+        if delta_l < delta_h
+            if alpha1==1 % no inner concavity
+                h = 1/gamma * (gamma * A(t) * alpha2/(w*(delta_h - delta_l)))^(1/(1-alpha2)) - n/gamma;
+                l = m - h;
+
+                if h < 0
+                    h = 0;
+                    l = m;
+                elseif h > m
+                    h = m;
+                    l = 0;
+                end
+            else % there is inner concavity
+                
+                if m < 1e-8
+                    h = 0;
+                    l = 0;
+                else
+                    
+                    int_h = linspace(1e-5, m, 100);
+                    trdoff = ((n^alpha1 + gamma * int_h.^alpha1).^(alpha2 - 1)).*(int_h.^(alpha1 - 1)) - ... 
+                               ((delta_h - delta_l) * w)/(A(t)*alpha1*alpha2*gamma);
+                    if trdoff(end) > 0
+                        h = m;
+                        l = 0;
+                    elseif trdoff(1) < 0
+                        h = 0;
+                        l = m;     
+                    else
+                        [min_trdoff, i_min_trdoff] = min(trdoff.^2);
+                        h = int_h(i_min_trdoff);
+                        l = m - h;     
+                        if ~isreal(h) || l < 0 || min_trdoff > 1e-5
+                            Trdoff = @(h) ((n^alpha1 + gamma * h^alpha1)^(alpha2 - 1))*(h^(alpha1 - 1)) - ... 
+                               ((delta_h - delta_l) * w)/(A(t)*alpha1*alpha2*gamma);
+                              h = fzero(Trdoff, [1e-5, m]); %Matlab's solver for h 
+                              l = m -  h; 
+                          if l < 0 || l > m
+                            error('out of bounds')
+                          end                         
+                        end
+                    end
+                    
+                end
+            end
+        else
+            l = 0;
+            h = m;
+        end
+        
+        if n < 0
+            n = 0;
+        end
+        
+        % if q is internalized, computes it
+        if flag_internalize_q==1
+            q(t) = Pi_q * (a_t + n_tilde_t + m_tilde_t);
+        end
+        
+        % production in t
+        production_t = A(t) * (n^alpha1 + gamma * h^alpha1)^alpha2;
+        
+        % profit in t
+        profit_t = production_t - ...
+            delta_n(t) * w * n - delta_h * w * h - ...
+            delta_l * w * l - delta_s * w * z_t;
+        
+        % discounted profit
+        disc_profits = disc_profits + beta^(t-1) * profit_t;
+        
+        % planner's objective function
+        if planner_obj_n==1 || planner_obj_n==2
+            planner_obj_t = production_t - delta_d * (d_t - d_tm1);
+        elseif planner_obj_n==3
+            planner_obj_t = production_t; 
+        end
+        
+        disc_planner_obj = disc_planner_obj + beta^(t-1) * planner_obj_t;
+        
+        % stores some output variables if asked to
+        if nargout>1
+            % calculates p_t
+            a_n_t = (nn_t + nm_t) * (a_t/(N - r_s_t - z_t - d_t));
+            p_t = Pi_pq * q(t) + Pi_pn * (a_n_t + n_tilde_t);
+            p_t = min(p_t, 1);
+            
+            % stores variables in vectors
+            n_vec(t) = n;
+            h_vec(t) = h;
+            l_vec(t) = l;
+            a_vec(t) = a_t;
+            s_vec(t) = z_t;
+            r_s_vec(t) = r_s_t;
+            r_a_vec(t) = r_a_t; 
+            d_vec(t) = d_t;
+            nn_vec(t) = nn_t;
+            mn_vec(t) = mn_t;
+            nm_vec(t) = nm_t;
+            mm_vec(t) = mm_t;
+            nr_vec(t) = nr_t;
+            mr_vec(t) = mr_t;
+            profit_vec(t) = profit_t;
+            infected_workplace_vec(t) = infected_workplace;
+            p_vec(t) = p_t;
+            n_tilde_vec(t) = n_tilde_t;
+            m_tilde_vec(t) = m_tilde_t;
+            production_vec(t) = production_t;
+        end
+        
+    end
+    
+    % profit in the no-disease world
+    delta_n_nd = 1;
+    production_nd = A(T) * (n_nd^alpha1 + gamma * m_nd^alpha1)^alpha2;
+    disc_profit_nd =  production_nd - delta_n_nd * w * n_nd - delta_h * w * m_nd;
+    disc_profit_nd = disc_profit_nd/(1-beta);
+    
+    % same for planner
+    disc_planner_obj_nd = production_nd/(1-beta);
+    
+    % variables in t=T+1
+    n_Tp1 = (N - d_t - z_t) * n_nd/N;
+    m_Tp1 = (N - d_t - z_t) * m_nd/N;
+    production_Tp1 = A(T) * (n_Tp1^alpha1 + gamma * m_Tp1^alpha1)^alpha2;
+    profit_Tp1 = production_Tp1 - delta_n_nd * w * n_Tp1 - delta_h * w * m_Tp1;
+    
+    % adds future discounted profits after T
+    disc_profits = disc_profits + beta^T * profit_Tp1 / (1-beta);
+    
+    % same for planner
+    if planner_obj_n==1
+        disc_planner_obj = disc_planner_obj + beta^T * production_Tp1 / (1-beta);
+    elseif planner_obj_n==3
+        disc_planner_obj = disc_planner_obj + beta^T * production_Tp1 / (1-beta) - ...
+            delta_d * d_t;
+    end
+    
+%     disc_planner_obj = disc_planner_obj + beta^T * profit_Tp1 / (1-beta);
+    
+%     disc_planner_obj = disc_planner_obj + beta^T * production_Tp1 / (1-beta);
+%     disc_planner_obj = delta_d * (- d_t) + (1 - delta_d) * disc_planner_obj;
+    
+    if flag_planner==0
+        out = (disc_profits/disc_profit_nd - 1)*1e+7;
+    else
+        out = (disc_planner_obj/disc_planner_obj_nd - 1)*1e+7;
+    end
+    
+    % we are looking for a minimum
+    out = - out;
+    
+    % stores output variables in a structdigits(digitsOld) if asked to
+    if nargout>1
+        % some variables
+        incubation = n_tilde_vec + m_tilde_vec;
+        exposed = [incubation(2:T); incubation(T)];
+        susceptible = N - s_vec - a_vec - r_s_vec - r_a_vec - d_vec - ...
+            incubation - exposed;
+        
+        % saves variables to a struct
+        outstruct = struct( ...
+            'q', q, ...
+            'n', n_vec, ...
+            'h', h_vec, ...
+            'l', l_vec, ...
+            'a', a_vec, ...
+            's', s_vec, ...
+            'r_s', r_s_vec, ...
+            'r_a', r_a_vec, ...
+            'd', d_vec, ...
+            'nn', nn_vec, ...
+            'mn', mn_vec, ...
+            'nm', nm_vec, ...
+            'mm', mm_vec, ...
+            'nr', nr_vec, ...
+            'mr', mr_vec, ...
+            'profit', profit_vec, ...
+            'infected_workplace', infected_workplace_vec, ...
+            'p', p_vec, ...
+            'n_tilde', n_tilde_vec, ...
+            'm_tilde', m_tilde_vec, ...
+            'production', production_vec, ...
+            'incubation', incubation, ...
+            'exposed', exposed, ...
+            'susceptible', susceptible, ...
+            'A', A, ...
+            'disc_profits', disc_profits, ...
+            'solution', in);
+    end
+    
+end
+
+
+
+
+
+
diff --git a/models/BJ_22/equilibrium.m b/models/BJ_22/equilibrium.m
new file mode 100644
index 0000000..0853931
--- /dev/null
+++ b/models/BJ_22/equilibrium.m
@@ -0,0 +1,139 @@
+
+
+% initial guess for q
+if flag_fix_q==1
+    q = q_fixed;
+    q_m1_fix = q_m1_fixed;
+else
+    q = zeros(T, 1);
+    q(1:20) = linspace(0, 0.04, 20);
+    q(21:40) = linspace(0.04, 0, 20);
+
+    q_m1 = Pi_q * initial_infected;
+end
+
+% guess for choices
+if flag_use_guess==0
+    if alpha1==1
+        guess = ones(T, 2) * N;
+    else
+        guess = ones(T, 3) * (1 - h_over_n_plus_h) * N;
+    end
+end
+
+% bounds for choices
+if alpha1==1
+    lb = zeros(T, 2);
+    ub = ones(T, 2) * N;
+else
+    lb = zeros(T, 3);
+    ub = ones(T, 3) * N;
+end
+
+% general equilibrium tolerance
+tol_ge = 1e-6; %1e-7;
+norm = tol_ge + 1;
+
+% general equilibrium loop
+while norm > tol_ge
+    
+    % gets firm's choices
+    if flag_fixed_choices==1
+        % fixed choices
+        solution = guess;
+    else
+        % solves problem of firm
+%         solution = fmincon(@compute_profits, guess, [], [], [], [], lb, ub, [], opts_fmincon);
+        
+%         opts = optimset('MaxFunEvals', 1e10, 'MaxIter', 1e10);
+%         solution = fminsearch(@compute_profits, guess, opts);
+
+max_iter = intmax;
+tol = 1e-10;
+
+opts_fmincon = optimoptions(@fmincon, ...
+          'MaxIter', 1e10, ...
+          'MaxFunEvals', 1e10, ...
+          'FunctionTolerance', tol, ...
+          'OptimalityTolerance', tol, ...
+          'StepTolerance', tol, ...
+          'ConstraintTolerance', tol, ...
+          'TolCon', tol, ...
+          'Display', 'off');
+
+%         opts = optimset('MaxFunEvals', 1e10, 'MaxIter', 1e10);
+%         solution = fmincon(@compute_profits, guess, [],[],[],[],lb,ub,[], opts);
+        problem = createOptimProblem('fmincon',...
+                'objective',@(x)compute_profits(x),...
+                'x0', guess, 'lb', lb, 'ub', ub, ...
+                'options', opts_fmincon);
+        solution = fmincon(problem);
+%                 problem = createOptimProblem('fmincon',...
+%                 'objective',@(x)compute_profits(x),...
+%                 'x0', solution, 'lb', lb, 'ub', ub, ...
+%                 'options', opts_fmincon);
+%         gs = GlobalSearch('Display','off');
+%         ms = MultiStart('Display','off');
+%         solution = run(gs, problem);
+%         solution = run(ms,problem, 200);
+%         solution = patternsearch(@compute_profits, guess, [],[],[],[],lb,ub,[], opts);
+    end
+    
+    % use solution of this iteration as guess for next iteration
+    guess = solution;
+    
+    % gets n. of infectious agents to compute new q
+    [~, outstruct] = compute_profits(solution);
+    
+    % new q
+    if flag_fix_q==1
+        q_m1_new = q_m1_fixed;
+        q_new = q_fixed;
+    else
+        infectious = outstruct.a + outstruct.n_tilde + outstruct.m_tilde;
+        q_m1_new = Pi_q * initial_infected;
+        q_new = Pi_q * infectious;
+    end
+    
+    % defines A if it is endogenous
+    if flag_endog_A_1==1
+        s = outstruct.s;
+        s_max = max(s);
+        A_new = 1 + s * (A_min - 1)/s_max;
+    elseif flag_endog_A_2==1
+        s = outstruct.s;
+        s_max = max(s);
+        A_new = 1 + slope_A * s;
+    else
+        A_new = A;
+    end
+    
+    % defines delta_n if it is endogenous
+    if flag_endog_delta_n==1
+        s = outstruct.s;
+        s_max = max(s);
+        delta_n_new =  s/s_max * delta_n_max + (s_max - s)/s_max * 1;
+    else
+        delta_n_new = delta_n;
+    end
+    
+    % checks convergence of equilibrium variables
+    V = [q_m1; q; A; delta_n];
+    V_new = [q_m1_new; q_new; A_new; delta_n_new];
+    
+    norm = max(abs(V - V_new));
+    
+    % updates equilibrium variables
+    q_m1 = q_m1_new;
+    q = q_new;
+    A = A_new;
+    delta_n = delta_n_new;
+    
+    % display
+    fprintf('norm: %.6e\n', norm)
+    
+end
+
+
+
+
diff --git a/models/BJ_22/main.m b/models/BJ_22/main.m
new file mode 100644
index 0000000..c78f7c9
--- /dev/null
+++ b/models/BJ_22/main.m
@@ -0,0 +1,452 @@
+
+% clc
+% clear
+% close all
+
+% tic
+
+global N T q q_m1 rho_r_s rho_r_a w delta_n delta_h delta_s delta_l A ...
+    beta alpha1 alpha2 gamma rho_d initial_infected varphi Pi_pq Pi_pn ...
+    h_over_n_plus_h n_nd m_nd flag_internalize_q Pi_q flag_planner ...
+    delta_d planner_obj_n
+
+%-------------------------------------------------------------------------%
+%                                 Moments                                 %
+%-------------------------------------------------------------------------%
+
+R0_data = 2.5; % R0 for COVID
+h_over_n_plus_h = 1 - 0.9823667407; % h/(n+h) in the no-disease case
+
+%-------------------------------------------------------------------------%
+%                          Technical parameters                           %
+%-------------------------------------------------------------------------%
+
+% T = 120; % time horizon
+T = 52 * 3; % time horizon
+% T = 52 * 1.5; % time horizon
+
+% fmincon parameters
+max_iter = intmax;
+tol = 1e-10;
+
+% opts_fmincon = optimoptions('fmincon', ...
+%           'Algorithm', 'active-set', ...  
+%           'MaxIter', max_iter, ...
+%           'MaxFunEvals', max_iter, ...
+%           'FunctionTolerance', tol, ...
+%           'OptimalityTolerance', tol, ...
+%           'StepTolerance', tol, ...
+%           'ConstraintTolerance', tol, ...
+%           'TolCon', tol, ...
+%           'Display', 'None');
+
+%-------------------------------------------------------------------------%
+%                               Parameters                                %
+%-------------------------------------------------------------------------%
+
+alpha1 = 1; % production function concavity (inside)
+alpha2 = 0.7; % production function concavity (outside)
+gamma = 0.957; % teleworking productivity
+
+N = 1; % nr. of workers
+%initial_infected = 0.001; % initial number of infected workers
+helper=load('inf_ini.mat');
+initial_infected=helper.helper;
+beta = 0.96^(1/52); % time preference
+delta_n = ones(T, 1); % relative wage rate
+delta_h = 1; % relative wage rate
+delta_l = 1; % relative wage rate
+delta_s = 1; % relative wage rate
+rho_r_s = 0.5; % probability of recovery (symptomatic)
+rho_r_a = 0.5; % probability of recovery (asymptomatic)
+rho_d = 0.0025;
+varphi = 1/2; % fraction of symptomatic
+A = ones(T, 1); % tfp
+% A = [linspace(1, 0.9, T/25)'; linspace(0.9, 0.9, 2*T/25)'; linspace(0.9, 1, 9*T/25)'; ones(T * 15/25, 1)];
+% delta_n = [linspace(1, 5, T/25)'; linspace(5, 5, 2*T/25)'; linspace(5, 1, 9*T/25)'; ones(T * 13/25, 1)];
+
+Pi_pq = 1; % infection rate parameter (related to intercept of p)
+Pi_pn = 0.7;%0.75; %1/N * 2/3; % infection rate parameter (slope) (0.45/1.2) 2/3
+
+T_a = 1/rho_r_a;
+T_i = 1 + (1-varphi) * T_a;
+
+% Pi_q = (R0_data/((1-initial_infected) * T_i) - Pi_pn)/Pi_pq; % Pi_q to fit R0_data
+% fprintf('Pi_q: %.8f\n', Pi_q)
+
+Pi_q = 0.55; %1/4;
+
+Pi2 = Pi_q;
+
+% R0 = (1-initial_infected) * T_i * (Pi_q * Pi_pq + Pi_pn);
+R0 = T_i * (Pi_q * Pi_pq + Pi_pn);
+fprintf('R0: %.8f\n', R0)
+
+% qq = Pi_q * initial_infected;
+% pp = Pi_pq * qq + Pi_pn * initial_infected;
+% R0 = T_i * pp/initial_infected % transmissions per initially infected
+
+%-------------------------------------------------------------------------%
+%                Sets some parameters by matching moments                 %
+%-------------------------------------------------------------------------%
+
+if alpha1==1 % h is zero in the no-disease case
+    
+    h_over_n_plus_h = 0;
+    
+    w = A(T)*alpha2/(delta_n(T) * N^(1-alpha2)); % w s.t. optimal demand is N
+    
+    n_nd = N;
+    m_nd = 0;
+    
+    fprintf('w: %.8f\n', w)
+    
+else % there is inada for h
+    
+    % finds w such that n + h = N in no-COVID world
+    n_nd = (1-h_over_n_plus_h) * N;
+    m_nd = h_over_n_plus_h * N;
+    
+    w = A(T) * (n_nd^alpha1 + gamma * m_nd^alpha1)^(alpha2-1) * alpha1 * alpha2 * n_nd^(alpha1-1);
+    
+    fprintf('w: %.8f\n', w)
+
+end
+
+%-------------------------------------------------------------------------%
+%                               Some flags                                %
+%-------------------------------------------------------------------------%
+
+flag_fixed_choices = 0; % if 1, uses fixed choices for the firm
+
+flag_fix_q = 0; % if 1, uses q_fixed and q_m1_fixed for q and q_m1
+
+flag_endog_A_1 = 0; % if 1, A is a function of the number of sick workers
+                    % uses variable A_min
+
+A_min = 0.83;
+
+flag_endog_A_2 = 0; % if 1, A is a function of the number of sick workers
+                    % uses slope_A
+                    
+slope_A = 0;
+
+flag_endog_delta_n = 0; % if 1, delta_n is a function of the n. of sick
+                        % workers. Uses variable_delta_n_max
+
+delta_n_max = 1.01;
+
+flag_internalize_q = 0; % if 1, firm internalizes its impact on q
+
+flag_planner = 0; % if 1, uses planner's objective function
+
+delta_d = 0;
+
+flag_use_guess = 0;
+
+planner_obj_n = 1;
+
+%-------------------------------------------------------------------------%
+%                              Runs something                             %
+%-------------------------------------------------------------------------%
+
+% variable that says what the code will do
+do = 1;
+
+if do==1 % finds general equilibrium
+    
+    % calculates general equilibrium
+    equilibrium
+    
+    % starts a cell
+    C = cell(1, 1);
+    
+    % saves output variables in this cell
+    [~, C{1}] = compute_profits(solution);
+    
+    % makes figures
+ %   i_fig = 1;
+%     filename_fig_suffix = 'benchmark'; % define this variable if you want to save figure
+%  figures_July9
+    
+elseif do==2 % simulates several scenarios
+    
+    policies
+    
+elseif do==3 % some figures
+    
+    % starts a cell
+    C = cell(5, 1);
+    
+    % equil 1: benchmark
+    equilibrium
+    [~, C{1}] = compute_profits(solution);
+    
+    % equil. 2: fixed choices
+    flag_fixed_choices = 1;
+    equilibrium
+    [~, C{2}] = compute_profits(solution);
+    flag_fixed_choices = 0;
+    
+    % equil. 3: delta_l = 0
+    delta_l_bm = delta_l;
+    delta_l = 0;
+    equilibrium
+    [~, C{3}] = compute_profits(solution);
+    delta_l = delta_l_bm;
+    
+    % equil. 4: A shock
+    A_min = A_min;
+    flag_endog_A_1 = 1;
+    equilibrium
+    [~, C{4}] = compute_profits(solution);
+    
+    % equil. 5: large A shock
+    A_min = 0.5;
+    flag_endog_A_1 = 1;
+    equilibrium
+    [~, C{5}] = compute_profits(solution);
+    
+    % makes figures
+    i_vec = [1, 3, 4, 5];
+    i_overlay_vec = [2, 1, 1, 1];
+    fig_desc_vec = {'benchmark', 'delta_l', 'tfp', 'tfplarge'};
+    
+    for ii = 1:length(i_vec)
+        for i_overlay = 0:1
+            if ii==1
+                leg_out = {'Benchmark', 'Fixed choices'};
+            elseif ii==4
+                leg_out = {'Large changes in $A$', 'Benchmark'};
+            else
+                leg_out = {'Counterfactual', 'Benchmark'};
+            end
+            
+            if i_overlay==0
+                i_fig = i_vec(ii);
+                leg_out = '';
+            else
+                i_fig = [i_vec(ii), i_overlay_vec(ii)];
+            end
+            
+            filename_fig_suffix = [fig_desc_vec{ii}, '_overlay', num2str(i_overlay)];
+            figures_July13
+            close all
+        end
+    end
+    
+elseif do==4 % TFP shock with benchmark q
+    
+    % starts a cell
+    C = cell(3, 1);
+    
+    % minimum A
+    A_min = 0.833;
+    flag_endog_A_1 = 0;
+    
+    % calculates general equilibrium
+    equilibrium
+    
+    % saves benchmark q
+    q_fixed = q;
+    q_m1_fixed = q_m1;
+    
+    % saves output variables
+    [~, C{1}] = compute_profits(solution);
+    
+    % TFP shock
+    s = C{1}.s;
+    s_max = max(s);
+    A = 1 + s * (A_min - 1)/s_max;
+    
+    % fixes benchmark q
+    flag_fix_q = 1;
+    
+    % finds choices of firm
+    equilibrium
+    
+    % saves output variables
+    [~, C{2}] = compute_profits(solution);
+    
+    % sets variables for equilibrium with TFP shock
+    flag_fix_q = 0;
+    flag_endog_A_1 = 1;
+    
+    % calculates general equilibrium
+    equilibrium
+    
+    % saves output variables
+    [~, C{3}] = compute_profits(solution);
+    
+    % adds some variables to structs
+    for i = 1:3
+        C{i}.incubation = C{i}.n_tilde + C{i}.m_tilde;
+        C{i}.exposed = [C{i}.incubation(2:T); C{i}.incubation(T)];
+        C{i}.susceptible = N - C{i}.s - C{i}.a - C{i}.r_s - C{i}.r_a - C{i}.d - C{i}.incubation - C{i}.exposed;
+        C{i}.available = N - C{i}.d - C{i}.s;
+    end
+    
+    % figures 1
+    i_fig = [2, 1];
+    leg_out = {'Counterfactual', 'Benchmark'};
+    t0_fig = 1;
+    t1_fig = 60;
+    figures_July13
+    
+    % figures 2
+    i_fig = [3, 2];
+    leg_out = {'General equilibrium', 'Partial equilibrium'};
+    t0_fig = 1;
+    t1_fig = 60;
+    figures_July13
+    
+    % figures 3
+    i_fig = [3, 1];
+    leg_out = {'TFP shock (g.e.)', 'Benchmark'};
+    t0_fig = 1;
+    t1_fig = 60;
+    figures_July13
+    
+elseif do==5 % policies
+    
+    policies
+    
+elseif do==6 % policies with endogenous A
+    
+    flag_endog_A_1 = 1;
+    A_min = 0.8315;
+    
+    policies
+    
+elseif do==7 % planner's problem (and the rest)
+    
+    planner_obj_n = 1;
+    
+    % vector with delta values
+    delta_d_vec = [0, 0, 0, 1e2, 1e3, 1e10];
+%     delta_d_vec = [0, 0, 0, 100];
+%     delta_d_vec = [701, 702, 703, 704];
+    
+    % starts a cell
+    C = cell(length(delta_d_vec), 1);
+    
+    % runs planner's problem for many deltas
+    for i_delta_d_vec = 1:length(delta_d_vec)
+        
+        if i_delta_d_vec==1 % benchmark
+            flag_internalize_q = 0;
+            flag_planner = 0;
+        elseif i_delta_d_vec==2 % benchmark, internalize q
+            flag_internalize_q = 1;
+            flag_planner = 0;
+        else
+            flag_internalize_q = 1;
+            flag_planner = 1;
+            
+            delta_d = delta_d_vec(i_delta_d_vec);
+        end
+        
+        % solves planner's problem
+        equilibrium
+        
+        % saves output variables in a cell
+        [~, C{i_delta_d_vec}] = compute_profits(solution);
+
+        % adds delta_d to C
+        C{i_delta_d_vec}.delta_d = delta_d;
+        
+        % use this solution as guess for next simulation problem
+%         flag_use_guess = 1;
+    end
+
+    filename = ['mat/planner_obj_n_1_T_', num2str(T), '.mat'];
+
+    save(filename)
+    
+elseif do==8 % different gamma and delta_n
+    
+
+    flag_internalize_q = 0;
+    flag_planner = 0;
+    
+    % vector with delta values
+    gamma = linspace(0.70, 0.97, 20);
+    delta_n_vec = linspace(1, 1.30, 20);    
+    param_vec = combvec(gamma, delta_n_vec);
+    
+    
+    % starts a cell
+    C = cell(length(param_vec), 1);
+    
+    % runs planner's problem for many deltas
+    for i_param_vec = 1:length(param_vec)
+        
+        gamma = param_vec(1, i_param_vec);
+        delta_n = param_vec(2, i_param_vec).* ones(T, 1);        
+        % solves planner's problem
+        equilibrium
+        
+        % saves output variables in a cell
+        [~, C{i_param_vec}] = compute_profits(solution);
+
+        % adds delta_d to C
+        C{i_param_vec}.delta_h = delta_h;
+        C{i_param_vec}.delta_n = delta_n;        
+        
+    end    
+
+%     mesh(linspace(0.70, 0.97, 10), delta_n_vec, reshape(X1(14,1:end), [C,C]))
+
+end
+
+
+
+
+
+ 
+% toc
+
+%-------------------------------------------------------------------------%
+%                                Functions                                %
+%-------------------------------------------------------------------------%
+
+% this function evaluates the FOC for choices of the firm in the no-disease
+% case
+function out = n_h_nd(in)
+    
+    global A alpha1 alpha2 gamma delta_h delta_n w T
+
+    % "in" is a guess for h
+
+    % transforms in to positive
+    h = in^2;
+    
+    % gets n
+    n = (delta_h/(gamma * delta_n))^(1/(1-alpha1)) * h;
+    
+    % evaluates FOC
+    out = A(T) * alpha2 * (n^alpha1 + gamma * h^alpha1)^(alpha2 - 1) * ...
+        gamma * alpha1 * h^(alpha1 - 1) - delta_h * w;
+
+end
+
+% this functions finds parameters such that optimal demand of firm in the
+% no-disease case fits h/(n+h) = n_over_n_plus_h and n + h = N
+function out = calib_nd(in)
+    
+    global gamma w delta_n delta_h N h_over_n_plus_h alpha1
+
+    % defines parameters
+    gamma = 0.5 * (1 + sin(in(1)));
+    w = in(2)^2;
+    
+    % finds optimal choices of n and h
+    h = fzero(@n_h_nd, 1);
+    h = h^2;
+    n = (delta_h/(gamma * delta_n))^(1/(1-alpha1)) * h;
+    
+    % evaluates conditions to be fitted
+    out(1) = h/(n+h) - h_over_n_plus_h;
+    out(2) = n + h - N;
+    
+end
\ No newline at end of file
diff --git a/models/BJ_22/modelmasterscript.m b/models/BJ_22/modelmasterscript.m
new file mode 100644
index 0000000..9f6aa1b
--- /dev/null
+++ b/models/BJ_22/modelmasterscript.m
@@ -0,0 +1,26 @@
+main;
+
+Consumption = NaN(156,1);
+Labour = C{1}.n+C{1}.h;
+Output = C{1}.production;
+Susceptibles = C{1}.susceptible; 
+Infected = C{1}.exposed + C{1}.incubation + C{1}.s + C{1}.a;
+Recovered = C{1}.r_a + C{1}.r_s;
+Deaths = C{1}.d;
+Interest = NaN(length(Consumption),1);
+Inflation = NaN(length(Consumption),1);
+Investment = NaN(length(Consumption),1);
+
+Consumption_ss = 0;
+Labour_ss = C{1}.n(1)+C{1}.h(1);
+Output_ss = C{1}.production(1);
+Susceptibles_ss= 0;
+Infected_ss= 0;
+Recovered_ss= 0;
+Deaths_ss= 0;
+Interest_ss = 0;
+Inflation_ss = 0;
+Investment_ss = 0;
+
+save('simulated_results.mat','Consumption','Labour','Output','Deaths','Susceptibles','Infected','Recovered','Interest','Inflation','Investment');
+save('simulated_results_ss.mat','Consumption_ss','Labour_ss','Output_ss','Deaths_ss','Susceptibles_ss','Infected_ss','Recovered_ss','Interest_ss','Inflation_ss','Investment_ss');
diff --git a/models/BJ_22/simulated_results.mat b/models/BJ_22/simulated_results.mat
new file mode 100644
index 0000000..d5f7961
Binary files /dev/null and b/models/BJ_22/simulated_results.mat differ
diff --git a/models/BJ_22/simulated_results_ss.mat b/models/BJ_22/simulated_results_ss.mat
new file mode 100644
index 0000000..1c69993
Binary files /dev/null and b/models/BJ_22/simulated_results_ss.mat differ