diff --git a/code/calibration/__init__.py b/code/calibration/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/code/calibration/estimation_parameters.py b/code/calibration/estimation_parameters.py
index c195f8f..a2f596a 100644
--- a/code/calibration/estimation_parameters.py
+++ b/code/calibration/estimation_parameters.py
@@ -74,73 +74,7 @@
 )
 
 # Age-varying discount factors over the lifecycle, lifted from Cagetti (2003)
-DiscFac_timevary = [
-    1.064914,
-    1.057997,
-    1.051422,
-    1.045179,
-    1.039259,
-    1.033653,
-    1.028352,
-    1.023348,
-    1.018632,
-    1.014198,
-    1.010037,
-    1.006143,
-    1.002509,
-    0.9991282,
-    0.9959943,
-    0.9931012,
-    0.9904431,
-    0.9880143,
-    0.9858095,
-    0.9838233,
-    0.9820506,
-    0.9804866,
-    0.9791264,
-    0.9779656,
-    0.9769995,
-    0.9762239,
-    0.9756346,
-    0.9752274,
-    0.9749984,
-    0.9749437,
-    0.9750595,
-    0.9753422,
-    0.9757881,
-    0.9763936,
-    0.9771553,
-    0.9780698,
-    0.9791338,
-    0.9803439,
-    0.981697,
-    0.8287214,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-    0.9902111,
-]
+DiscFac_timevary = np.genfromtxt("code/data/Cagetti2003.csv")
 
 # Survival probabilities over the lifecycle
 liv_prb = parse_ssa_life_table(
@@ -150,15 +84,10 @@
 # Age groups for the estimation: calculate average wealth-to-permanent income ratio
 # for consumers within each of these age groups, compare actual to simulated data
 empirical_cohort_age_groups = [
-    [26, 27, 28, 29, 30],
-    [31, 32, 33, 34, 35],
-    [36, 37, 38, 39, 40],
-    [41, 42, 43, 44, 45],
-    [46, 47, 48, 49, 50],
-    [51, 52, 53, 54, 55],
-    [56, 57, 58, 59, 60],
+    [age for age in range(start, start + 5)] for start in range(26, 61, 5)
 ]
 
+
 # Three point discrete distribution of initial w
 initial_wealth_income_ratio_vals = np.array([0.17, 0.5, 0.83])
 # Equiprobable discrete distribution of initial w
diff --git a/code/calibration/setup_scf_data.py b/code/calibration/setup_scf_data.py
index 9d29484..fc071e1 100644
--- a/code/calibration/setup_scf_data.py
+++ b/code/calibration/setup_scf_data.py
@@ -2,44 +2,20 @@
 Sets up the SCF data for use in the SolvingMicroDSOPs estimation.
 """
 
-import os
+from code.calibration.estimation_parameters import (
+    empirical_cohort_age_groups,
+    initial_age,
+)
 
 import numpy as np  # Numerical Python
 import pandas as pd
 
-from calibration.estimation_parameters import empirical_cohort_age_groups, initial_age
 
-# Find pathname to this file:
-my_file_path = os.path.dirname(os.path.abspath(__file__))
-
-# Pathnames to the other files:
-# Relative directory for primitive parameter files
-calibration_dir = os.path.join(my_file_path, "../calibration/")
-# Relative directory for primitive parameter files
-tables_dir = os.path.join(my_file_path, "../tables/")
-# Relative directory for primitive parameter files
-figures_dir = os.path.join(my_file_path, "../figures/")
-# Relative directory for primitive parameter files
-code_dir = os.path.join(my_file_path, "../code/")
-
-# Need to rely on the manual insertion of pathnames to all files in do_all.py
-# NOTE sys.path.insert(0, os.path.abspath(tables_dir)), etc. may need to be
-# copied from do_all.py to here
-
-# Import files first:
-
-
-# The following libraries are part of the standard python distribution
-
-# Set the path to the empirical data:
-# os.path.abspath('./')   #'./'
-scf_data_path = data_location = os.path.dirname(os.path.abspath(__file__))
-
-data = pd.read_csv(scf_data_path + "/data/SCFdata.csv")
+scf_data = pd.read_csv("code/data/SCFdata.csv")
 
 # Keep only observations with normal incomes > 0,
 # otherwise wealth/income is not well defined
-data = data[data.norminc > 0.0]
+scf_data = scf_data[scf_data.norminc > 0.0]
 
 # Initialize empty lists for the data
 w_to_y_data = []  # Ratio of wealth to permanent income
@@ -49,7 +25,7 @@
 # Only extract the data required from the SCF dataset
 search_ages = [ages[-1] for ages in empirical_cohort_age_groups]
 for idx, age in enumerate(search_ages):
-    age_data = data[data.age_group.str.contains(f"{age}]")]
+    age_data = scf_data[scf_data.age_group.str.contains(f"{age}]")]
     w_to_y_data.append(age_data.wealth_income_ratio.values)
     empirical_weights.append(age_data.weight.values)
     # create a group id 1-7 for every observation
diff --git a/code/data/Cagetti2003.csv b/code/data/Cagetti2003.csv
new file mode 100644
index 0000000..e769b01
--- /dev/null
+++ b/code/data/Cagetti2003.csv
@@ -0,0 +1,65 @@
+1.064914
+1.057997
+1.051422
+1.045179
+1.039259
+1.033653
+1.028352
+1.023348
+1.018632
+1.014198
+1.010037
+1.006143
+1.002509
+0.9991282
+0.9959943
+0.9931012
+0.9904431
+0.9880143
+0.9858095
+0.9838233
+0.9820506
+0.9804866
+0.9791264
+0.9779656
+0.9769995
+0.9762239
+0.9756346
+0.9752274
+0.9749984
+0.9749437
+0.9750595
+0.9753422
+0.9757881
+0.9763936
+0.9771553
+0.9780698
+0.9791338
+0.9803439
+0.981697
+0.8287214
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
+0.9902111
diff --git a/code/calibration/data/SCFdata.csv b/code/data/SCFdata.csv
similarity index 100%
rename from code/calibration/data/SCFdata.csv
rename to code/data/SCFdata.csv
diff --git a/do_all.py b/code/do_all.py
similarity index 81%
rename from do_all.py
rename to code/do_all.py
index cf5c41b..8f9f75b 100644
--- a/do_all.py
+++ b/code/do_all.py
@@ -46,43 +46,14 @@
 still run.
 """
 
-import os
-import sys
-
-from calibration.options import (
+from code.calibration.options import (
+    all_replications,
+    high_resource,
     low_resource,
     medium_resource,
-    high_resource,
-    all_replications,
 )
 from code.estimation import estimate
 
-# Find pathname to this file:
-my_file_path = os.path.dirname(os.path.abspath(__file__))
-
-# Pathnames to the other files:
-# Relative directory for primitive parameter files
-calibration_dir = os.path.join(my_file_path, "calibration")
-# Relative directory for primitive parameter files
-tables_dir = os.path.join(my_file_path, "tables")
-# Relative directory for primitive parameter files
-figures_dir = os.path.join(my_file_path, "figures")
-# Relative directory for primitive parameter files
-code_dir = os.path.join(my_file_path, "code")
-
-# manually add the pathnames to the various files directly to the beginning
-# of the Python path. This will be needed for all files that will run in
-# lower directories.
-sys.path.insert(0, calibration_dir)
-sys.path.insert(0, tables_dir)
-sys.path.insert(0, figures_dir)
-sys.path.insert(0, code_dir)
-sys.path.insert(0, my_file_path)
-
-
-# Manual import needed, should draw from first instance at start of Python
-# PATH added above:
-
 
 # Ask the user which replication to run, and run it:
 def run_replication():
diff --git a/do_figs.py b/code/do_figs.py
similarity index 100%
rename from do_figs.py
rename to code/do_figs.py
diff --git a/code/estimation.py b/code/estimation.py
index 6a3740a..d03fa8c 100644
--- a/code/estimation.py
+++ b/code/estimation.py
@@ -9,9 +9,10 @@
 income as defined in ConsIndShockModel.
 """
 
+# Parameters for the consumer type and the estimation
+import code.calibration.estimation_parameters as parameters
+import code.calibration.setup_scf_data as data  # SCF 2004 data on household wealth
 import csv
-import os
-import sys
 from time import time  # Timing utility
 
 import matplotlib.pyplot as plt
@@ -35,29 +36,16 @@
 from HARK.estimation import bootstrap_sample_from_data, minimize_nelder_mead
 from scipy.optimize import approx_fprime
 
-# Parameters for the consumer type and the estimation
-import calibration.estimation_parameters as parameters
-import calibration.setup_scf_data as data  # SCF 2004 data on household wealth
-
-# Find pathname to this file:
-my_file_path = os.path.dirname(os.path.abspath(__file__))
 
 # Pathnames to the other files:
 # Relative directory for primitive parameter files
-calibration_dir = os.path.join(my_file_path, "../calibration/")
+calibration_dir = "code/calibration/"
 # Relative directory for primitive parameter files
-tables_dir = os.path.join(my_file_path, "../tables/")
+tables_dir = "content/tables/"
 # Relative directory for primitive parameter files
-figures_dir = os.path.join(my_file_path, "../figures/")
+figures_dir = "content/figures/"
 # Relative directory for primitive parameter files
-code_dir = os.path.join(my_file_path, "../code/")
-
-# Add the calibration folder to the path
-sys.path.insert(0, os.path.abspath(calibration_dir))
-
-# Need to rely on the manual insertion of pathnames to all files in do_all.py
-# NOTE sys.path.insert(0, os.path.abspath(tables_dir)), etc. may need to be
-# copied from do_all.py to here
+code_dir = "code/"
 
 
 # Set booleans to determine which tasks should be done
@@ -513,9 +501,10 @@ def smmObjectiveFxnReduced(parameters_to_estimate):
         )
 
         # Create the simple estimate table
-        estimate_results_file = os.path.join(
-            tables_dir, estimation_agent + "_estimate_results.csv"
+        estimate_results_file = (
+            tables_dir + "/" + estimation_agent + "_estimate_results.csv"
         )
+
         with open(estimate_results_file, "wt") as f:
             writer = csv.writer(f)
             writer.writerow(["DiscFacAdj", "CRRA"])
@@ -559,9 +548,10 @@ def smmObjectiveFxnReduced(parameters_to_estimate):
         )
 
         # Create the simple bootstrap table
-        bootstrap_results_file = os.path.join(
-            tables_dir, estimation_agent + "_bootstrap_results.csv"
+        bootstrap_results_file = (
+            tables_dir + "/" + estimation_agent + "_bootstrap_results.csv"
         )
+
         with open(bootstrap_results_file, "wt") as f:
             writer = csv.writer(f)
             writer.writerow(
@@ -630,9 +620,9 @@ def simulate_moments_reduced(x):
         axs[1].set_ylabel("Sensitivity")
         axs[1].set_xlabel("Median W/Y Ratio")
 
-        plt.savefig(os.path.join(figures_dir, estimation_agent + "Sensitivity.pdf"))
-        plt.savefig(os.path.join(figures_dir, estimation_agent + "Sensitivity.png"))
-        plt.savefig(os.path.join(figures_dir, estimation_agent + "Sensitivity.svg"))
+        plt.savefig(figures_dir + "/" + estimation_agent + "Sensitivity.pdf")
+        plt.savefig(figures_dir + "/" + estimation_agent + "Sensitivity.png")
+        plt.savefig(figures_dir + "/" + estimation_agent + "Sensitivity.svg")
 
         plt.show()
 
@@ -675,9 +665,9 @@ def simulate_moments_reduced(x):
         pylab.plot(model_estimate[1], model_estimate[0], "*r", ms=15)
         pylab.xlabel(r"coefficient of relative risk aversion $\rho$", fontsize=14)
         pylab.ylabel(r"discount factor adjustment $\beth$", fontsize=14)
-        pylab.savefig(os.path.join(figures_dir, estimation_agent + "SMMcontour.pdf"))
-        pylab.savefig(os.path.join(figures_dir, estimation_agent + "SMMcontour.png"))
-        pylab.savefig(os.path.join(figures_dir, estimation_agent + "SMMcontour.svg"))
+        pylab.savefig(figures_dir + "/" + estimation_agent + "SMMcontour.pdf")
+        pylab.savefig(figures_dir + "/" + estimation_agent + "SMMcontour.png")
+        pylab.savefig(figures_dir + "/" + estimation_agent + "SMMcontour.svg")
         pylab.show()
 
 
@@ -871,14 +861,14 @@ def simulate_moments_reduced(x):
             f"Time to execute all: {time_to_contour / 60:.2f} min, {time_to_contour:.2f} sec"
         )
 
-    fig_sensitivity.savefig(os.path.join(figures_dir, "AllSensitivity.pdf"))
-    fig_sensitivity.savefig(os.path.join(figures_dir, "AllSensitivity.png"))
-    fig_sensitivity.savefig(os.path.join(figures_dir, "AllSensitivity.svg"))
+    fig_sensitivity.savefig(figures_dir + "/AllSensitivity.pdf")
+    fig_sensitivity.savefig(figures_dir + "/AllSensitivity.png")
+    fig_sensitivity.savefig(figures_dir + "/AllSensitivity.svg")
 
     # Save and show the plot
-    fig_contour.savefig(os.path.join(figures_dir, "AllSMMcontour.pdf"))
-    fig_contour.savefig(os.path.join(figures_dir, "AllSMMcontour.png"))
-    fig_contour.savefig(os.path.join(figures_dir, "AllSMMcontour.svg"))
+    fig_contour.savefig(figures_dir + "/AllSMMcontour.pdf")
+    fig_contour.savefig(figures_dir + "/AllSMMcontour.png")
+    fig_contour.savefig(figures_dir + "/AllSMMcontour.svg")
 
 
 if __name__ == "__main__":
diff --git a/code/tests.py b/code/tests.py
index f4d371e..8e82378 100644
--- a/code/tests.py
+++ b/code/tests.py
@@ -1,4 +1,4 @@
-from calibration.options import low_resource, medium_resource
+from code.calibration.options import low_resource, medium_resource
 from code.estimation import estimate
 
 
diff --git a/content/paper/math.ipynb b/content/paper/math.ipynb
index 0223302..ff2eb48 100644
--- a/content/paper/math.ipynb
+++ b/content/paper/math.ipynb
@@ -7,7 +7,7 @@
    "outputs": [],
    "source": [
     "from sympy import *\n",
-    "from sympy.plotting import plot\n"
+    "from sympy.plotting import plot"
    ]
   },
   {
@@ -16,8 +16,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "c, a , x= symbols('c a x')\n",
-    "rho, delta = symbols('rho delta')"
+    "c, a, x = symbols(\"c a x\")\n",
+    "rho, delta = symbols(\"rho delta\")"
    ]
   },
   {
@@ -40,7 +40,7 @@
     }
    ],
    "source": [
-    "u = (c**(1-delta) * a**delta)**(1-rho)/(1-rho)\n",
+    "u = (c ** (1 - delta) * a**delta) ** (1 - rho) / (1 - rho)\n",
     "u"
    ]
   },
@@ -65,7 +65,7 @@
    ],
    "source": [
     "uc = simplify(u.diff(c))\n",
-    "uc\n"
+    "uc"
    ]
   },
   {
@@ -136,7 +136,9 @@
     }
    ],
    "source": [
-    "g = ((f*c*a/(-a * delta + a - c * delta))**(1/(1-rho))/a**delta)**(1/(1-delta))\n",
+    "g = ((f * c * a / (-a * delta + a - c * delta)) ** (1 / (1 - rho)) / a**delta) ** (\n",
+    "    1 / (1 - delta)\n",
+    ")\n",
     "simplify(g)"
    ]
   },
@@ -160,7 +162,7 @@
     }
    ],
    "source": [
-    "g_1 = g.subs({a:1000, delta:0.5, rho:2})\n",
+    "g_1 = g.subs({a: 1000, delta: 0.5, rho: 2})\n",
     "simplify(g_1)"
    ]
   },
@@ -207,7 +209,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "util = x**(1-rho)/(1-rho)"
+    "util = x ** (1 - rho) / (1 - rho)"
    ]
   },
   {
@@ -230,7 +232,7 @@
     }
    ],
    "source": [
-    "simplify(simplify(util.diff(x)**(-1/rho)))"
+    "simplify(simplify(util.diff(x) ** (-1 / rho)))"
    ]
   },
   {