diff --git a/GCN Files/RBC_two_household.gcn b/GCN Files/RBC_two_household.gcn
new file mode 100644
index 0000000..51bac29
--- /dev/null
+++ b/GCN Files/RBC_two_household.gcn
@@ -0,0 +1,177 @@
+assumptions
+{
+ positive
+ {
+ Y[], K[], C_NR[], C_R[],
+ w[], r[], mc_L[],
+ L[], L_NR[], L_R[],
+ TFP[],
+ alpha, alpha_L, beta, sigma_C, sigma_L, delta;
+ };
+};
+
+tryreduce
+{
+ U_NR[], U_R[], TC[];
+};
+
+block STEADY_STATE
+{
+ definitions
+ {
+ f1[ss] = (r[ss] / (r[ss] - alpha * delta));
+ f2[ss] = (alpha_L * (1 - alpha_L) * (1 - alpha)) ^ (-sigma_L / sigma_C);
+ f3[ss] = alpha_L ^ (sigma_L / sigma_C) +
+ (1 - alpha_L) ^ (sigma_L / sigma_C);
+
+ };
+ identities
+ {
+ TFP[ss] = 1.0;
+ shock_beta_R[ss] = 1.0;
+
+ r[ss] = 1 / beta - (1 - delta);
+ w[ss] = (1 - alpha) * alpha_L ^ alpha_L * (1 - alpha_L) ^ (1 - alpha_L) *
+ (r[ss] / alpha) ^ (alpha / (alpha - 1));
+ mc_L[ss] = w[ss] / alpha_L ^ alpha_L / (1 - alpha_L) ^ (1 - alpha_L);
+ Y[ss] = (f1[ss] * f2[ss] * f3[ss] * w[ss] ^ ((1 + sigma_L) / sigma_C)) ^
+ (sigma_C / (sigma_L + sigma_C));
+
+ L[ss] = (1 - alpha) * Y[ss] / mc_L[ss];
+
+ L_R[ss] = alpha_L * L[ss] * mc_L[ss] / w[ss];
+ L_NR[ss] = (1 - alpha_L) * L[ss] * mc_L[ss] / w[ss];
+
+ C_R[ss] = w[ss] ^ (1/sigma_C) * L_R[ss] ^ (-sigma_L / sigma_C);
+ C_NR[ss] = w[ss] ^ (1/sigma_C) * L_NR[ss] ^ (-sigma_L / sigma_C);
+
+ K[ss] = alpha * Y[ss] / r[ss];
+ I[ss] = delta * K[ss];
+
+ lambda_R[ss] = C_R[ss] ^ -sigma_C;
+ q[ss] = lambda_R[ss];
+ lambda_NR[ss] = C_NR[ss] ^ -sigma_C;
+
+ };
+};
+
+block RICARDIAN_HOUSEHOLD
+{
+ definitions
+ {
+ u_R[] = shock_beta_R[] * (C_R[] ^ (1 - sigma_C) / (1 - sigma_C) -
+ L_R[] ^ (1 + sigma_L) / (1 + sigma_L));
+ };
+
+ controls
+ {
+ C_R[], L_R[], I[], K[];
+ };
+
+ objective
+ {
+ U_R[] = u_R[] + beta * E[][U_R[1]];
+ };
+
+ constraints
+ {
+ @exclude
+ C_R[] + I[] = r[] * K[-1] + w[] * L_R[] : lambda_R[];
+
+ K[] = (1 - delta) * K[-1] + I[]: q[];
+ };
+
+ identities
+ {
+ log(shock_beta_R[]) = rho_beta_R * log(shock_beta_R[-1]) + epsilon_beta_R[];
+ };
+
+ shocks
+ {
+ epsilon_beta_R[];
+ };
+
+ calibration
+ {
+ beta = 0.99;
+ delta = 0.02;
+ sigma_C = 1.5;
+ sigma_L = 2.0;
+ rho_beta_R = 0.95;
+ };
+};
+
+block NON_RICARDIAN_HOUSEHOLD
+{
+ definitions
+ {
+ u_NR[] = (C_NR[] ^ (1 - sigma_C) / (1 - sigma_C) -
+ L_NR[] ^ (1 + sigma_L) / (1 + sigma_L));
+ };
+
+ controls
+ {
+ C_NR[], L_NR[];
+ };
+
+ objective
+ {
+ U_NR[] = u_NR[] + beta * E[][U_NR[1]];
+ };
+
+ constraints
+ {
+ @exclude
+ C_NR[] = w[] * L_NR[]: lambda_NR[];
+ };
+};
+
+
+block FIRM
+{
+ controls
+ {
+ K[-1], L[], L_R[], L_NR[];
+ };
+
+ objective
+ {
+ TC[] = -(r[] * K[-1] + w[] * L_R[] + w[] * L_NR[]);
+ };
+
+ constraints
+ {
+ L[] = L_R[] ^ alpha_L * L_NR[] ^ (1 - alpha_L) : mc_L[];
+ Y[] = TFP[] * K[-1] ^ alpha * L[] ^ (1 - alpha) : mc[];
+ };
+
+ identities
+ {
+ # Perfect competition
+ mc[] = 1;
+
+ # Exogenous technology process
+ log(TFP[]) = rho_TFP * log(TFP[-1]) + epsilon_TFP[];
+ };
+
+ shocks
+ {
+ epsilon_TFP[];
+ };
+
+ calibration
+ {
+ alpha = 0.35;
+ alpha_L = 0.5;
+
+ rho_TFP = 0.95;
+ };
+};
+
+block EQULIBRIUM
+{
+ identities
+ {
+ Y[] = C_R[] + C_NR[] + I[];
+ };
+};
diff --git a/GCN Files/RBC_two_household_additive.txt b/GCN Files/RBC_two_household_additive.txt
new file mode 100644
index 0000000..636aefb
--- /dev/null
+++ b/GCN Files/RBC_two_household_additive.txt
@@ -0,0 +1,183 @@
+assumptions
+{
+ positive
+ {
+ Y[], K[], C_NR[], C_R[],
+ w[], r[], mc_L[],
+ L[], L_NR[], L_R[],
+ TFP[],
+ alpha, alpha_L, beta, sigma_C, sigma_L, delta, omega;
+ };
+};
+
+tryreduce
+{
+ U_NR[], U_R[], TC[];
+};
+
+# block STEADY_STATE
+# {
+ # definitions
+ # {
+ # f1[ss] = (r[ss] / (r[ss] - alpha * delta));
+ # f2[ss] = (alpha_L * (1 - alpha_L) * (1 - alpha)) ^ (-sigma_L / sigma_C);
+ # f3[ss] = alpha_L ^ (sigma_L / sigma_C) +
+ # (1 - alpha_L) ^ (sigma_L / sigma_C);
+
+ # };
+ # identities
+ # {
+ # TFP[ss] = 1.0;
+ # shock_beta_R[ss] = 1.0;
+
+ # r[ss] = 1 / beta - (1 - delta);
+ # w[ss] = (1 - alpha) * alpha_L ^ alpha_L * (1 - alpha_L) ^ (1 - alpha_L) *
+ # (r[ss] / alpha) ^ (alpha / (alpha - 1));
+ # mc_L[ss] = w[ss] / alpha_L ^ alpha_L / (1 - alpha_L) ^ (1 - alpha_L);
+ # Y[ss] = (f1[ss] * f2[ss] * f3[ss] * w[ss] ^ ((1 + sigma_L) / sigma_C)) ^
+ # (sigma_C / (sigma_L + sigma_C));
+
+ # L[ss] = (1 - alpha) * Y[ss] / mc_L[ss];
+
+ # L_R[ss] = alpha_L * L[ss] * mc_L[ss] / w[ss];
+ # L_NR[ss] = (1 - alpha_L) * L[ss] * mc_L[ss] / w[ss];
+
+ # C_R[ss] = w[ss] ^ (1/sigma_C) * L_R[ss] ^ (-sigma_L / sigma_C);
+ # C_NR[ss] = w[ss] ^ (1/sigma_C) * L_NR[ss] ^ (-sigma_L / sigma_C);
+
+ # K[ss] = alpha * Y[ss] / r[ss];
+ # I[ss] = delta * K[ss];
+
+ # lambda_R[ss] = C_R[ss] ^ -sigma_C;
+ # q[ss] = lambda_R[ss];
+ # lambda_NR[ss] = C_NR[ss] ^ -sigma_C;
+
+ # };
+# };
+
+block RICARDIAN_HOUSEHOLD
+{
+ definitions
+ {
+ u_R[] = shock_beta_R[] * (C_R[] ^ (1 - sigma_C) / (1 - sigma_C) -
+ L_R[] ^ (1 + sigma_L) / (1 + sigma_L));
+
+ Phi_I[] = (1 - gamma_I / 2 * (I[] / I[-1] - 1) ^ 2);
+ };
+
+ controls
+ {
+ C_R[], L_R[], I[], K[];
+ };
+
+ objective
+ {
+ U_R[] = u_R[] + beta * E[][U_R[1]];
+ };
+
+ constraints
+ {
+ @exclude
+ C_R[] + I[] = r[] * K[-1] + w[] * L_R[] : lambda_R[];
+
+ K[] = (1 - delta) * K[-1] + Phi_I[] * I[]: q[];
+ };
+
+ identities
+ {
+ log(shock_beta_R[]) = rho_beta_R * log(shock_beta_R[-1]) + epsilon_beta_R[];
+ };
+
+ shocks
+ {
+ epsilon_beta_R[];
+ };
+
+ calibration
+ {
+ beta = 0.99;
+ delta = 0.02;
+ sigma_C = 1.5;
+ sigma_L = 2.0;
+ rho_beta_R = 0.95;
+ gamma_I = 5;
+ };
+};
+
+block NON_RICARDIAN_HOUSEHOLD
+{
+ definitions
+ {
+ u_NR[] = (C_NR[] ^ (1 - sigma_C) / (1 - sigma_C) -
+ L_NR[] ^ (1 + sigma_L) / (1 + sigma_L));
+ };
+
+ controls
+ {
+ C_NR[], L_NR[];
+ };
+
+ objective
+ {
+ U_NR[] = u_NR[] + beta * E[][U_NR[1]];
+ };
+
+ constraints
+ {
+ C_NR[] = w[] * L_NR[]: lambda_NR[];
+ };
+};
+
+
+block FIRM
+{
+ controls
+ {
+ K[-1], L[];
+ };
+
+ objective
+ {
+ TC[] = -(r[] * K[-1] + w[] * L[]);
+ };
+
+ constraints
+ {
+ Y[] = TFP[] * K[-1] ^ alpha * L[] ^ (1 - alpha) : mc[];
+ };
+
+ identities
+ {
+ # Perfect competition
+ mc[] = 1;
+
+ # Exogenous technology process
+ log(TFP[]) = rho_TFP * log(TFP[-1]) + epsilon_TFP[];
+ };
+
+ shocks
+ {
+ epsilon_TFP[];
+ };
+
+ calibration
+ {
+ alpha = 0.35;
+ rho_TFP = 0.95;
+ };
+};
+
+block EQULIBRIUM
+{
+ identities
+ {
+ L[] = omega * L_R[] + (1 - omega) * L_NR[];
+ C[] = omega * C_R[] + (1 - omega) * C_NR[];
+ Y[] = C[] + I[];
+ };
+
+ calibration
+ {
+ omega = 0.5;
+ };
+};
diff --git a/GCN Files/skilled_unskilled_rbc.gcn b/GCN Files/skilled_unskilled_rbc.gcn
new file mode 100644
index 0000000..797f3ce
--- /dev/null
+++ b/GCN Files/skilled_unskilled_rbc.gcn
@@ -0,0 +1,134 @@
+block STEADY_STATE
+{
+ identities
+ {
+ A[ss] = 1.0;
+ Div[ss] = 0.0;
+ r_u[ss] = 1 / beta_u - (1 - delta_u);
+ r_s[ss] = 1 / beta_s - (1 - delta_s);
+ };
+};
+
+block SKILLED_HOUSEHOLD
+{
+ definitions
+ {
+ u_s[] = log(C_s[]) - Theta_s * L_s[];
+ };
+
+ objective
+ {
+ U_s[] = u_s[] + beta_s * E[][U_s[1]];
+ };
+
+ controls
+ {
+ C_s[], L_s[], K_s[], I_s[];
+ };
+
+ constraints
+ {
+ C_s[] + I_s[] = w_s[] * L_s[] + r_s[] * K_s[-1] + s * Div[]: lambda_s[];
+ K_s[] = (1 - delta_s) * K_s[-1] + I_s[];
+ };
+
+ calibration
+ {
+ beta_s = 0.99;
+ delta_s = 0.035;
+ Theta_s = 1;
+ s = 0.5; # Share of dividend that the skilled household gets (could be alpha_L ?)
+ };
+};
+
+block UNSKILLED_HOUSEHOLD
+{
+ definitions
+ {
+ u_u[] = log(C_u[]) - Theta_u * L_u[];
+ };
+
+ objective
+ {
+ U_u[] = u_u[] + beta_u * E[][U_u[1]];
+ };
+
+ controls
+ {
+ C_u[], L_u[], K_u[], I_u[];
+ };
+
+ constraints
+ {
+ C_u[] + I_u[] = w_u[] * L_u[] + r_u[] * K_u[-1] + (1 - s) * Div[]: lambda_u[];
+ K_u[] = (1 - delta_u) * K_u[-1] + I_u[];
+ };
+
+ calibration
+ {
+ beta_u = 0.99;
+ delta_u = 0.035;
+ Theta_u = 1;
+ };
+};
+
+
+block FIRM
+{
+ objective
+ {
+ TC[] = -(r_u[] * K_u[] + r_s[] * K_s[] + w_u[] * L_u[] + w_s[] * L_s[]);
+ };
+
+ controls
+ {
+ K_u[-1], K_s[-1], L_u[], L_s[], K[], L[];
+ };
+
+ constraints
+ {
+ # Bundle labor -- skilled/unskilled are imperfect substitutes
+ L[] = (alpha_L ^ (1 / psi_L) * L_u[] ^ ((psi_L - 1) / psi_L) +
+ (1 - alpha_L) ^ (1 / psi_L) * L_s[] ^ ((psi_L - 1) / psi_L)) ^
+ (psi_L / (psi_L - 1));
+
+ # Bundle capital -- perfect substitutes
+ K[] = K_u[-1] ^ alpha_K * K_s[-1] ^ (1 - alpha_K);
+
+ # Production function
+ Y[] = A[] * K[] ^ alpha * L[] ^ (1 - alpha) : P[];
+ };
+
+ identities
+ {
+ # Perfect competition
+ P[] = 1;
+ Div[] = Y[] * P[] + TC[];
+ };
+
+ calibration
+ {
+ alpha_L = 0.5; # share of unskilled labor in economy
+ alpha_K = 0.5; # share of capital stock owned by unskilled household
+ psi_L = 3.0; # Elasticity of substitution btwn skilled & unskilled, psi_L -> oo implies perfect substitutes
+ alpha = 0.66; # Share of capital in production
+ };
+};
+
+block TECHNOLOGY
+{
+ identities
+ {
+ log(A[]) = rho * log(A[-1]) + epsilon_A[];
+ };
+
+ shocks
+ {
+ epsilon_A[];
+ };
+
+ calibration
+ {
+ rho = 0.95;
+ };
+};
diff --git a/examples/Multiple Households.ipynb b/examples/Multiple Households.ipynb
new file mode 100644
index 0000000..7d48c05
--- /dev/null
+++ b/examples/Multiple Households.ipynb
@@ -0,0 +1,420 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 181,
+ "id": "a3c80085",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import sympy as sp\n",
+ "import gEconpy as ge\n",
+ "import gEconpy.plotting as gp"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 258,
+ "id": "844c4ca8",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Model Building Complete.\n",
+ "Found:\n",
+ "\t16 equations\n",
+ "\t16 variables\n",
+ "\tThe following variables were eliminated at user request:\n",
+ "\t\tTC_t,U_NR_t,U_R_t\n",
+ "\tThe following \"variables\" were defined as constants and have been substituted away:\n",
+ "\t\tmc_t\n",
+ "\t2 stochastic shocks\n",
+ "\t\t 0 / 2 has a defined prior. \n",
+ "\t9 parameters\n",
+ "\t\t 0 / 9 has a defined prior. \n",
+ "\t0 parameters to calibrate.\n",
+ "Model appears well defined and ready to proceed to solving.\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "mod = ge.model_from_gcn(\"../../../gEconpy/GCN Files/RBC_two_household_additive.gcn\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 262,
+ "id": "cf73612b",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "976cc10634da4f9eac50342e9245f146",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "Output()"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n"
+ ],
+ "text/plain": []
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Steady state found\n",
+ "--------------------------------------------------------------------------------\n",
+ "Optimizer message The solution converged.\n",
+ "Sum of squared residuals 1.4317932293331621e-21\n",
+ "Maximum absoluate error 2.69482214321215e-11\n",
+ "Gradient L2-norm at solution 1.7978379435964637e-10\n",
+ "Max abs gradient at solution 9.149458968238378e-11\n"
+ ]
+ }
+ ],
+ "source": [
+ "ss_res, success = mod.steady_state(how=\"root\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 271,
+ "id": "8e61081b",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "Model solution has 4 eigenvalues greater than one in modulus and 4 forward-looking variables. \n",
+ "Blanchard-Kahn condition is satisfied.\n"
+ ]
+ },
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Modulus | \n",
+ " Real | \n",
+ " Imaginary | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 1.653468e-35 | \n",
+ " -1.653468e-35 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 1.653028e-34 | \n",
+ " -1.653028e-34 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " 1.686148e-34 | \n",
+ " 1.686148e-34 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " 7.255700e-23 | \n",
+ " -7.255700e-23 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " 7.609220e-23 | \n",
+ " -7.609220e-23 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 5 | \n",
+ " 9.155729e-22 | \n",
+ " 9.155729e-22 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 6 | \n",
+ " 1.014493e-21 | \n",
+ " -1.014493e-21 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 7 | \n",
+ " 3.309565e-21 | \n",
+ " 3.309565e-21 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 8 | \n",
+ " 4.360073e-21 | \n",
+ " 4.360073e-21 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 9 | \n",
+ " 1.318637e-19 | \n",
+ " 1.318637e-19 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 10 | \n",
+ " 1.439463e-17 | \n",
+ " -1.439463e-17 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 11 | \n",
+ " 2.998660e-17 | \n",
+ " 2.998660e-17 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 12 | \n",
+ " 7.283610e-01 | \n",
+ " -7.283610e-01 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 13 | \n",
+ " 9.500000e-01 | \n",
+ " -9.500000e-01 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 14 | \n",
+ " 9.500000e-01 | \n",
+ " -9.500000e-01 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 15 | \n",
+ " 9.735574e-01 | \n",
+ " -9.735574e-01 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 16 | \n",
+ " 1.027303e+00 | \n",
+ " -1.027303e+00 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 17 | \n",
+ " 1.400627e+00 | \n",
+ " -1.400627e+00 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 18 | \n",
+ " 4.756147e+06 | \n",
+ " 4.756147e+06 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ " | 19 | \n",
+ " 4.922745e+07 | \n",
+ " 4.922745e+07 | \n",
+ " 0.0 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "