diff --git a/solutions/stand_alone_programs/bisection2.py b/solutions/stand_alone_programs/bisection2.py
deleted file mode 100644
index 958dd5404..000000000
--- a/solutions/stand_alone_programs/bisection2.py
+++ /dev/null
@@ -1,26 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: bisection2.py
-Authors: John Stachurski, Thomas J. Sargent
-LastModified: 11/08/2013
-
-"""
-
-def bisect(f, a, b, tol=10e-5):
-    """
-    Implements the bisection root finding algorithm, assuming that f is a
-    real-valued function on [a, b] satisfying f(a) < 0 < f(b).
-    """
-    lower, upper = a, b
-    if upper - lower < tol:
-        return 0.5 * (upper + lower)
-    else:
-        middle = 0.5 * (upper + lower)
-        print('Current mid point = {}'.format(middle))
-        if f(middle) > 0:   # Implies root is between lower and middle
-            bisect(f, lower, middle)
-        else:               # Implies root is between middle and upper
-            bisect(f, middle, upper)
-
-
-
diff --git a/solutions/stand_alone_programs/discrete_rv0.py b/solutions/stand_alone_programs/discrete_rv0.py
deleted file mode 100644
index 537bbf971..000000000
--- a/solutions/stand_alone_programs/discrete_rv0.py
+++ /dev/null
@@ -1,33 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: discrete_rv0.py
-Authors: John Stachurski and Thomas Sargent
-LastModified: 11/08/2013
-
-"""
-
-from numpy import cumsum
-from numpy.random import uniform
-
-class discreteRV:
-    """
-    Generates an array of draws from a discrete random variable with vector of
-    probabilities given by q.  
-    """
-
-    def __init__(self, q):
-        """
-        The argument q is a NumPy array, or array like, nonnegative and sums
-        to 1
-        """
-        self.q = q
-        self.Q = cumsum(q)
-
-    def draw(self, k=1):
-        """
-        Returns k draws from q. For each such draw, the value i is returned
-        with probability q[i].
-        """
-        return self.Q.searchsorted(uniform(0, 1, size=k)) 
-
-
diff --git a/solutions/stand_alone_programs/graph.txt b/solutions/stand_alone_programs/graph.txt
deleted file mode 100644
index 51dc5ba7e..000000000
--- a/solutions/stand_alone_programs/graph.txt
+++ /dev/null
@@ -1,100 +0,0 @@
-node0, node1 0.04, node8 11.11, node14 72.21
-node1, node46 1247.25, node6 20.59, node13 64.94
-node2, node66 54.18, node31 166.80, node45 1561.45
-node3, node20 133.65, node6 2.06, node11 42.43
-node4, node75 3706.67, node5 0.73, node7 1.02
-node5, node45 1382.97, node7 3.33, node11 34.54
-node6, node31 63.17, node9 0.72, node10 13.10
-node7, node50 478.14, node9 3.15, node10 5.85
-node8, node69 577.91, node11 7.45, node12 3.18
-node9, node70 2454.28, node13 4.42, node20 16.53
-node10, node89 5352.79, node12 1.87, node16 25.16
-node11, node94 4961.32, node18 37.55, node20 65.08
-node12, node84 3914.62, node24 34.32, node28 170.04
-node13, node60 2135.95, node38 236.33, node40 475.33
-node14, node67 1878.96, node16 2.70, node24 38.65
-node15, node91 3597.11, node17 1.01, node18 2.57
-node16, node36 392.92, node19 3.49, node38 278.71
-node17, node76 783.29, node22 24.78, node23 26.45
-node18, node91 3363.17, node23 16.23, node28 55.84
-node19, node26 20.09, node20 0.24, node28 70.54
-node20, node98 3523.33, node24 9.81, node33 145.80
-node21, node56 626.04, node28 36.65, node31 27.06
-node22, node72 1447.22, node39 136.32, node40 124.22
-node23, node52 336.73, node26 2.66, node33 22.37
-node24, node66 875.19, node26 1.80, node28 14.25
-node25, node70 1343.63, node32 36.58, node35 45.55
-node26, node47 135.78, node27 0.01, node42 122.00
-node27, node65 480.55, node35 48.10, node43 246.24
-node28, node82 2538.18, node34 21.79, node36 15.52
-node29, node64 635.52, node32 4.22, node33 12.61
-node30, node98 2616.03, node33 5.61, node35 13.95
-node31, node98 3350.98, node36 20.44, node44 125.88
-node32, node97 2613.92, node34 3.33, node35 1.46
-node33, node81 1854.73, node41 3.23, node47 111.54
-node34, node73 1075.38, node42 51.52, node48 129.45
-node35, node52 17.57, node41 2.09, node50 78.81
-node36, node71 1171.60, node54 101.08, node57 260.46
-node37, node75 269.97, node38 0.36, node46 80.49
-node38, node93 2767.85, node40 1.79, node42 8.78
-node39, node50 39.88, node40 0.95, node41 1.34
-node40, node75 548.68, node47 28.57, node54 53.46
-node41, node53 18.23, node46 0.28, node54 162.24
-node42, node59 141.86, node47 10.08, node72 437.49
-node43, node98 2984.83, node54 95.06, node60 116.23
-node44, node91 807.39, node46 1.56, node47 2.14
-node45, node58 79.93, node47 3.68, node49 15.51
-node46, node52 22.68, node57 27.50, node67 65.48
-node47, node50 2.82, node56 49.31, node61 172.64
-node48, node99 2564.12, node59 34.52, node60 66.44
-node49, node78 53.79, node50 0.51, node56 10.89
-node50, node85 251.76, node53 1.38, node55 20.10
-node51, node98 2110.67, node59 23.67, node60 73.79
-node52, node94 1471.80, node64 102.41, node66 123.03
-node53, node72 22.85, node56 4.33, node67 88.35
-node54, node88 967.59, node59 24.30, node73 238.61
-node55, node84 86.09, node57 2.13, node64 60.80
-node56, node76 197.03, node57 0.02, node61 11.06
-node57, node86 701.09, node58 0.46, node60 7.01
-node58, node83 556.70, node64 29.85, node65 34.32
-node59, node90 820.66, node60 0.72, node71 0.67
-node60, node76 48.03, node65 4.76, node67 1.63
-node61, node98 1057.59, node63 0.95, node64 4.88
-node62, node91 132.23, node64 2.94, node76 38.43
-node63, node66 4.43, node72 70.08, node75 56.34
-node64, node80 47.73, node65 0.30, node76 11.98
-node65, node94 594.93, node66 0.64, node73 33.23
-node66, node98 395.63, node68 2.66, node73 37.53
-node67, node82 153.53, node68 0.09, node70 0.98
-node68, node94 232.10, node70 3.35, node71 1.66
-node69, node99 247.80, node70 0.06, node73 8.99
-node70, node76 27.18, node72 1.50, node73 8.37
-node71, node89 104.50, node74 8.86, node91 284.64
-node72, node76 15.32, node84 102.77, node92 133.06
-node73, node83 52.22, node76 1.40, node90 243.00
-node74, node81 1.07, node76 0.52, node78 8.08
-node75, node92 68.53, node76 0.81, node77 1.19
-node76, node85 13.18, node77 0.45, node78 2.36
-node77, node80 8.94, node78 0.98, node86 64.32
-node78, node98 355.90, node81 2.59
-node79, node81 0.09, node85 1.45, node91 22.35
-node80, node92 121.87, node88 28.78, node98 264.34
-node81, node94 99.78, node89 39.52, node92 99.89
-node82, node91 47.44, node88 28.05, node93 11.99
-node83, node94 114.95, node86 8.75, node88 5.78
-node84, node89 19.14, node94 30.41, node98 121.05
-node85, node97 94.51, node87 2.66, node89 4.90
-node86, node97 85.09
-node87, node88 0.21, node91 11.14, node92 21.23
-node88, node93 1.31, node91 6.83, node98 6.12
-node89, node97 36.97, node99 82.12
-node90, node96 23.53, node94 10.47, node99 50.99
-node91, node97 22.17
-node92, node96 10.83, node97 11.24, node99 34.68
-node93, node94 0.19, node97 6.71, node99 32.77
-node94, node98 5.91, node96 2.03
-node95, node98 6.17, node99 0.27
-node96, node98 3.32, node97 0.43, node99 5.87
-node97, node98 0.30
-node98, node99 0.33
-node99, 
diff --git a/solutions/stand_alone_programs/linapprox.py b/solutions/stand_alone_programs/linapprox.py
deleted file mode 100644
index 5ffe4e79d..000000000
--- a/solutions/stand_alone_programs/linapprox.py
+++ /dev/null
@@ -1,44 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: linapprox.py
-Authors: John Stachurski, Thomas J. Sargent
-LastModified: 11/08/2013
-
-"""
-
-from __future__ import division  # Omit if using Python 3.x
-
-def linapprox(f, a, b, n, x):
-    """
-    Evaluates the piecewise linear interpolant of f at x on the interval 
-    [a, b], with n evenly spaced grid points.
-
-    Parameters 
-    ===========
-        f : function
-            The function to approximate
-
-        x, a, b : scalars (floats or integers) 
-            Evaluation point and endpoints, with a <= x <= b
-
-        n : integer
-            Number of grid points
-
-    Returns
-    =========
-        A float. The interpolant evaluated at x
-
-    """
-    length_of_interval = b - a
-    num_subintervals = n - 1
-    step = length_of_interval / num_subintervals  
-
-    # === find first grid point larger than x === #
-    point = a
-    while point <= x:
-        point += step
-
-    # === x must lie between the gridpoints (point - step) and point === #
-    u, v = point - step, point  
-
-    return f(u) + (x - u) * (f(v) - f(u)) / (v - u)
diff --git a/solutions/stand_alone_programs/solution_career_ex1.py b/solutions/stand_alone_programs/solution_career_ex1.py
deleted file mode 100644
index b6c1dd71b..000000000
--- a/solutions/stand_alone_programs/solution_career_ex1.py
+++ /dev/null
@@ -1,46 +0,0 @@
-""" 
-Simulate job / career paths and compute the waiting time to permanent job /
-career.  In reading the code, recall that optimal_policy[i, j] = policy at
-(theta_i, epsilon_j) = either 1, 2 or 3; meaning 'stay put', 'new job' and
-'new life'.
-"""
-import matplotlib.pyplot as plt
-import numpy as np
-from quantecon.discrete_rv import discreteRV
-from career import *
-from quantecon.compute_fp import compute_fixed_point
-
-wp = workerProblem()
-v_init = np.ones((wp.N, wp.N))*100
-v = compute_fixed_point(bellman, wp, v_init)
-optimal_policy = get_greedy(wp, v)
-F = discreteRV(wp.F_probs)
-G = discreteRV(wp.G_probs)
-
-def gen_path(T=20):
-    i = j = 0  
-    theta_index = []
-    epsilon_index = []
-    for t in range(T):
-        if optimal_policy[i, j] == 1:    # Stay put
-            pass
-        elif optimal_policy[i, j] == 2:  # New job
-            j = int(G.draw())
-        else:                            # New life
-            i, j  = int(F.draw()), int(G.draw())
-        theta_index.append(i)
-        epsilon_index.append(j)
-    return wp.theta[theta_index], wp.epsilon[epsilon_index]
-
-theta_path, epsilon_path = gen_path()
-fig = plt.figure()
-ax1 = plt.subplot(211)
-ax1.plot(epsilon_path, label='epsilon')
-ax1.plot(theta_path, label='theta')
-ax1.legend(loc='lower right')
-theta_path, epsilon_path = gen_path()
-ax2 = plt.subplot(212)
-ax2.plot(epsilon_path, label='epsilon')
-ax2.plot(theta_path, label='theta')
-ax2.legend(loc='lower right')
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_career_ex2.py b/solutions/stand_alone_programs/solution_career_ex2.py
deleted file mode 100644
index 6f4b6fa94..000000000
--- a/solutions/stand_alone_programs/solution_career_ex2.py
+++ /dev/null
@@ -1,32 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-from quantecon.discrete_rv import discreteRV
-from career import *
-from quantecon.compute_fp import compute_fixed_point
-
-wp = workerProblem()
-v_init = np.ones((wp.N, wp.N))*100
-v = compute_fixed_point(bellman, wp, v_init)
-optimal_policy = get_greedy(wp, v)
-F = discreteRV(wp.F_probs)
-G = discreteRV(wp.G_probs)
-
-def gen_first_passage_time():
-    t = 0
-    i = j = 0  
-    theta_index = []
-    epsilon_index = []
-    while 1:
-        if optimal_policy[i, j] == 1:    # Stay put
-            return t
-        elif optimal_policy[i, j] == 2:  # New job
-            j = int(G.draw())
-        else:                            # New life
-            i, j  = int(F.draw()), int(G.draw())
-        t += 1
-
-M = 25000 # Number of samples
-samples = np.empty(M)
-for i in range(M): 
-    samples[i] = gen_first_passage_time()
-print np.median(samples)
diff --git a/solutions/stand_alone_programs/solution_career_ex3.py b/solutions/stand_alone_programs/solution_career_ex3.py
deleted file mode 100644
index 34a7ea716..000000000
--- a/solutions/stand_alone_programs/solution_career_ex3.py
+++ /dev/null
@@ -1,22 +0,0 @@
-import matplotlib.pyplot as plt
-from matplotlib import cm
-from quantecon.career import *
-from compute_fp import compute_fixed_point
-
-wp = workerProblem()
-v_init = np.ones((wp.N, wp.N))*100
-v = compute_fixed_point(bellman, wp, v_init)
-optimal_policy = get_greedy(wp, v)
-
-fig = plt.figure(figsize=(6,6))
-ax = fig.add_subplot(111)
-tg, eg = np.meshgrid(wp.theta, wp.epsilon)
-lvls=(0.5, 1.5, 2.5, 3.5)
-ax.contourf(tg, eg, optimal_policy.T, levels=lvls, cmap=cm.winter, alpha=0.5)
-ax.contour(tg, eg, optimal_policy.T, colors='k', levels=lvls, linewidths=2)
-ax.set_xlabel('theta', fontsize=14)
-ax.set_ylabel('epsilon', fontsize=14)
-ax.text(1.8, 2.5, 'new life', fontsize=14)
-ax.text(4.5, 2.5, 'new job', fontsize=14, rotation='vertical')
-ax.text(4.0, 4.5, 'stay put', fontsize=14)
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_estspec_ex1.py b/solutions/stand_alone_programs/solution_estspec_ex1.py
deleted file mode 100644
index 16d4bd286..000000000
--- a/solutions/stand_alone_programs/solution_estspec_ex1.py
+++ /dev/null
@@ -1,33 +0,0 @@
-
-import numpy as np
-import matplotlib.pyplot as plt
-# from quantecon.linproc import linearProcess
-# from quantecon.estspec import periodogram
-import quantecon as qe
-
-
-## Data
-n = 400
-phi = 0.5
-theta = 0, -0.8
-lp = qe.linearProcess(phi, theta)
-X = lp.simulation(ts_length=n)
-
-fig, ax = plt.subplots(3, 1)
-
-for i, wl in enumerate((15, 55, 175)):  # window lengths
-    
-    x, y = qe.periodogram(X)
-    ax[i].plot(x, y, 'b-', lw=2, alpha=0.5, label='periodogram')
-
-    x_sd, y_sd = lp.spectral_density(two_pi=False, resolution=120)
-    ax[i].plot(x_sd, y_sd, 'r-', lw=2, alpha=0.8, label='spectral density')
-
-    x, y_smoothed = qe.periodogram(X, window='hamming', window_len=wl)
-    ax[i].plot(x, y_smoothed, 'k-', lw=2, label='smoothed periodogram')
-
-    ax[i].legend()
-    ax[i].set_title('window length = {}'.format(wl))
-
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_estspec_ex2.py b/solutions/stand_alone_programs/solution_estspec_ex2.py
deleted file mode 100644
index a025437e6..000000000
--- a/solutions/stand_alone_programs/solution_estspec_ex2.py
+++ /dev/null
@@ -1,26 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-import quantecon as qe
-
-lp = qe.linearProcess(-0.9)
-wl = 65
-
-
-fig, ax = plt.subplots(3, 1)
-
-for i in range(3):
-    X = lp.simulation(ts_length=150)
-    ax[i].set_xlim(0, np.pi)
-
-    x_sd, y_sd = lp.spectral_density(two_pi=False, resolution=180)
-    ax[i].semilogy(x_sd, y_sd, 'r-', lw=2, alpha=0.75, label='spectral density')
-
-    x, y_smoothed = qe.periodogram(X, window='hamming', window_len=wl)
-    ax[i].semilogy(x, y_smoothed, 'k-', lw=2, alpha=0.75, label='standard smoothed periodogram')
-
-    x, y_ar = qe.ar_periodogram(X, window='hamming', window_len=wl)
-    ax[i].semilogy(x, y_ar, 'b-', lw=2, alpha=0.75, label='AR smoothed periodogram')
-
-    ax[i].legend(loc='upper left')
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_ifp_ex1.py b/solutions/stand_alone_programs/solution_ifp_ex1.py
deleted file mode 100644
index 5ae46cc07..000000000
--- a/solutions/stand_alone_programs/solution_ifp_ex1.py
+++ /dev/null
@@ -1,29 +0,0 @@
-from matplotlib import pyplot as plt
-from quantecon.ifp import *
-
-m = consumerProblem()
-K = 80
-
-# Bellman iteration 
-V, c = initialize(m)
-print "Starting value function iteration"
-for i in range(K):
-    print "Current iterate = " + str(i)
-    V = bellman_operator(m, V)  
-c1 = bellman_operator(m, V, return_policy=True)  
-
-# Policy iteration 
-print "Starting policy function iteration"
-V, c2 = initialize(m)
-for i in range(K):
-    print "Current iterate = " + str(i)
-    c2 = coleman_operator(m, c2)
-
-fig, ax = plt.subplots()
-ax.plot(m.asset_grid, c1[:, 0], label='value function iteration')
-ax.plot(m.asset_grid, c2[:, 0], label='policy function iteration')
-ax.set_xlabel('asset level')
-ax.set_ylabel('consumption (low income)')
-ax.legend(loc='upper left')
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_ifp_ex2.py b/solutions/stand_alone_programs/solution_ifp_ex2.py
deleted file mode 100644
index 8ba590867..000000000
--- a/solutions/stand_alone_programs/solution_ifp_ex2.py
+++ /dev/null
@@ -1,19 +0,0 @@
-from compute_fp import compute_fixed_point
-from matplotlib import pyplot as plt
-import numpy as np
-from quantecon.ifp import coleman_operator, consumerProblem, initialize
-
-r_vals = np.linspace(0, 0.04, 4)  
-
-fig, ax = plt.subplots()
-for r_val in r_vals:
-    cp = consumerProblem(r=r_val)
-    v_init, c_init = initialize(cp)
-    c = compute_fixed_point(coleman_operator, cp, c_init)
-    ax.plot(cp.asset_grid, c[:, 0], label=r'$r = %.3f$' % r_val)
-
-ax.set_xlabel('asset level')
-ax.set_ylabel('consumption (low income)')
-ax.legend(loc='upper left')
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_ifp_ex3.py b/solutions/stand_alone_programs/solution_ifp_ex3.py
deleted file mode 100644
index c97d9839d..000000000
--- a/solutions/stand_alone_programs/solution_ifp_ex3.py
+++ /dev/null
@@ -1,33 +0,0 @@
-from matplotlib import pyplot as plt
-import numpy as np
-from quantecon.ifp import consumerProblem, coleman_operator, initialize
-from quantecon.compute_fp import compute_fixed_point
-from scipy import interp
-import mc_tools 
-
-def compute_asset_series(cp, T=500000):
-    """
-    Simulates a time series of length T for assets, given optimal savings
-    behavior.  Parameter cp is an instance of consumerProblem
-    """
-
-    Pi, z_vals, R = cp.Pi, cp.z_vals, cp.R  # Simplify names
-    v_init, c_init = initialize(cp)
-    c = compute_fixed_point(coleman_operator, cp, c_init)
-    cf = lambda a, i_z: interp(a, cp.asset_grid, c[:, i_z])
-    a = np.zeros(T+1)
-    z_seq = mc_tools.sample_path(Pi, sample_size=T)
-    for t in range(T):
-        i_z = z_seq[t]
-        a[t+1] = R * a[t] + z_vals[i_z] - cf(a[t], i_z)
-    return a
-
-if __name__ == '__main__':
-
-    cp = consumerProblem(r=0.03, grid_max=4)
-    a = compute_asset_series(cp)
-    fig, ax = plt.subplots()
-    ax.hist(a, bins=20, alpha=0.5, normed=True)
-    ax.set_xlabel('assets')
-    ax.set_xlim(-0.05, 0.75)
-    plt.show()
diff --git a/solutions/stand_alone_programs/solution_ifp_ex4.py b/solutions/stand_alone_programs/solution_ifp_ex4.py
deleted file mode 100644
index 12305e12e..000000000
--- a/solutions/stand_alone_programs/solution_ifp_ex4.py
+++ /dev/null
@@ -1,26 +0,0 @@
-from matplotlib import pyplot as plt
-import numpy as np
-from quantecon.compute_fp import compute_fixed_point
-from quantecon.ifp import coleman_operator, consumerProblem, initialize
-from solution_ifp_ex3 import compute_asset_series
-
-M = 25
-r_vals = np.linspace(0, 0.04, M)  
-fig, ax = plt.subplots()
-
-for b in (1, 3):
-    asset_mean = []
-    for r_val in r_vals:
-        cp = consumerProblem(r=r_val, b=b)
-        mean = np.mean(compute_asset_series(cp, T=250000))
-        asset_mean.append(mean)
-    ax.plot(asset_mean, r_vals, label=r'$b = %d$' % b)
-
-ax.set_yticks(np.arange(.0, 0.045, .01))
-ax.set_xticks(np.arange(-3, 2, 1))
-ax.set_xlabel('capital')
-ax.set_ylabel('interest rate')
-ax.grid(True)
-ax.legend(loc='upper left')
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_jv_ex1.py b/solutions/stand_alone_programs/solution_jv_ex1.py
deleted file mode 100644
index e1a690629..000000000
--- a/solutions/stand_alone_programs/solution_jv_ex1.py
+++ /dev/null
@@ -1,40 +0,0 @@
-import matplotlib.pyplot as plt
-import random
-from quantecon.jv import workerProblem, bellman_operator
-from quantecon.compute_fp import compute_fixed_point
-import numpy as np
-
-# Set up
-wp = workerProblem(grid_size=25)
-G, pi, F = wp.G, wp.pi, wp.F       # Simplify names
-
-v_init = wp.x_grid * 0.5
-V = compute_fixed_point(bellman_operator, wp, v_init, max_iter=40)
-s_policy, phi_policy = bellman_operator(wp, V, return_policies=True)
-
-# Turn the policy function arrays into actual functions
-s = lambda y: np.interp(y, wp.x_grid, s_policy)
-phi = lambda y: np.interp(y, wp.x_grid, phi_policy)
-
-def h(x, b, U):
-    return (1 - b) * G(x, phi(x)) + b * max(G(x, phi(x)), U)
-
-plot_grid_max, plot_grid_size = 1.2, 100
-plot_grid = np.linspace(0, plot_grid_max, plot_grid_size)
-fig, ax = plt.subplots()
-ax.set_xlim(0, plot_grid_max)
-ax.set_ylim(0, plot_grid_max)
-ticks = (0.25, 0.5, 0.75, 1.0)
-ax.set_xticks(ticks)
-ax.set_yticks(ticks)
-ax.set_xlabel(r'$x_t$', fontsize=16)
-ax.set_ylabel(r'$x_{t+1}$', fontsize=16, rotation='horizontal')
-
-ax.plot(plot_grid, plot_grid, 'k--')  # 45 degree line
-for x in plot_grid:
-    for i in range(50):
-        b = 1 if random.uniform(0, 1) < pi(s(x)) else 0
-        U = wp.F.rvs(1)
-        y = h(x, b, U)
-        ax.plot(x, y, 'go', alpha=0.25)
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_jv_ex2.py b/solutions/stand_alone_programs/solution_jv_ex2.py
deleted file mode 100644
index df42c2de7..000000000
--- a/solutions/stand_alone_programs/solution_jv_ex2.py
+++ /dev/null
@@ -1,16 +0,0 @@
-from matplotlib import pyplot as plt
-from quantecon.jv import workerProblem 
-import numpy as np
-
-# Set up
-wp = workerProblem(grid_size=25)
-
-def xbar(phi):
-    return (wp.A * phi**wp.alpha)**(1 / (1 - wp.alpha))
-
-phi_grid = np.linspace(0, 1, 100)
-fig, ax = plt.subplots()
-ax.set_xlabel(r'$\phi$', fontsize=16)
-ax.plot(phi_grid, [xbar(phi) * (1 - phi) for phi in phi_grid], 'b-', label=r'$w^*(\phi)$')
-ax.legend(loc='upper left')
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_kalman_ex1.py b/solutions/stand_alone_programs/solution_kalman_ex1.py
deleted file mode 100644
index e7a7ea2cf..000000000
--- a/solutions/stand_alone_programs/solution_kalman_ex1.py
+++ /dev/null
@@ -1,30 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.kalman import Kalman
-from scipy.stats import norm
-
-## Parameters
-theta = 10
-A, G, Q, R = 1, 1, 0, 1
-x_hat_0, Sigma_0 = 8, 1
-## Initialize Kalman filter
-kalman = Kalman(A, G, Q, R)
-kalman.set_state(x_hat_0, Sigma_0)
-
-N = 5
-fig, ax = plt.subplots()
-xgrid = np.linspace(theta - 5, theta + 2, 200)
-for i in range(N):
-    # Record the current predicted mean and variance, and plot their densities
-    m, v = kalman.current_x_hat, kalman.current_Sigma
-    m, v = float(m), float(v)
-    ax.plot(xgrid, norm.pdf(xgrid, loc=m, scale=np.sqrt(v)), label=r'$t=%d$' % i)
-    # Generate the noisy signal
-    y = theta + norm.rvs(size=1)
-    # Update the Kalman filter
-    kalman.update(y)
-
-ax.set_title(r'First %d densities when $\theta = %.1f$' % (N, theta)) 
-ax.legend(loc='upper left')
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_kalman_ex2.py b/solutions/stand_alone_programs/solution_kalman_ex2.py
deleted file mode 100644
index 8d0311ad9..000000000
--- a/solutions/stand_alone_programs/solution_kalman_ex2.py
+++ /dev/null
@@ -1,34 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.kalman import Kalman
-from scipy.stats import norm
-from scipy.integrate import quad
-
-## Parameters
-theta = 10
-A, G, Q, R = 1, 1, 0, 1
-x_hat_0, Sigma_0 = 8, 1
-epsilon = 0.1
-## Initialize Kalman filter
-kalman = Kalman(A, G, Q, R)
-kalman.set_state(x_hat_0, Sigma_0)
-
-T = 600
-z = np.empty(T)
-for t in range(T):
-    # Record the current predicted mean and variance, and plot their densities
-    m, v = kalman.current_x_hat, kalman.current_Sigma
-    m, v = float(m), float(v)
-    f = lambda x: norm.pdf(x, loc=m, scale=np.sqrt(v))
-    integral, error = quad(f, theta - epsilon, theta + epsilon)
-    z[t] = 1 - integral
-    # Generate the noisy signal and update the Kalman filter
-    kalman.update(theta + norm.rvs(size=1))
-
-fig, ax = plt.subplots()
-ax.set_ylim(0, 1)
-ax.set_xlim(0, T)
-ax.plot(range(T), z) 
-ax.fill_between(range(T), np.zeros(T), z, color="blue", alpha=0.2) 
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_kalman_ex3.py b/solutions/stand_alone_programs/solution_kalman_ex3.py
deleted file mode 100644
index 2c139c36f..000000000
--- a/solutions/stand_alone_programs/solution_kalman_ex3.py
+++ /dev/null
@@ -1,56 +0,0 @@
-from __future__ import print_function  # Remove for Python 3.x
-import numpy as np
-from numpy.random import multivariate_normal
-import matplotlib.pyplot as plt
-from scipy.linalg import eigvals
-from quantecon.kalman import Kalman
-
-# === Define A, Q, G, R === #
-G = np.eye(2)
-R = 0.5 * np.eye(2)
-A = [[0.5, 0.4], 
-     [0.6, 0.3]]
-Q = 0.3 * np.eye(2)
-
-# === Define the prior density === #
-Sigma = [[0.9, 0.3], 
-         [0.3, 0.9]]
-Sigma = np.array(Sigma)
-x_hat = np.array([8, 8])
-
-# === Initialize the Kalman filter === #
-kn = Kalman(A, G, Q, R)
-kn.set_state(x_hat, Sigma)
-
-# === Set the true initial value of the state === #
-x = np.zeros(2)
-
-# == Print eigenvalues of A == #
-print("Eigenvalues of A:")
-print(eigvals(A))
-
-# == Print stationary Sigma == #
-S, K = kn.stationary_values()
-print("Stationary prediction error variance:")
-print(S)
-
-# === Generate the plot === #
-T = 50
-e1 = np.empty(T)
-e2 = np.empty(T)
-for t in range(T):
-    # == Generate signal and update prediction == #
-    y = multivariate_normal(mean=np.dot(G, x), cov=R)
-    kn.update(y)
-    # == Update state and record error == #
-    Ax = np.dot(A, x)
-    x = multivariate_normal(mean=Ax, cov=Q)
-    e1[t] = np.sum((x - kn.current_x_hat)**2)
-    e2[t] = np.sum((x - Ax)**2)
-
-fig, ax = plt.subplots()
-ax.plot(range(T), e1, 'k-', lw=2, alpha=0.6, label='Kalman filter error') 
-ax.plot(range(T), e2, 'g-', lw=2, alpha=0.6, label='conditional expectation error') 
-ax.legend()
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_lln_ex1.py b/solutions/stand_alone_programs/solution_lln_ex1.py
deleted file mode 100644
index 533649260..000000000
--- a/solutions/stand_alone_programs/solution_lln_ex1.py
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-Illustrates the delta method, a consequence of the central limit theorem.
-"""
-import numpy as np
-from scipy.stats import uniform, norm
-import matplotlib.pyplot as plt
-from matplotlib import rc
-
-# == Specifying font, needs LaTeX integration == #
-rc('font',**{'family':'serif','serif':['Palatino']})
-rc('text', usetex=True)
-
-# == Set parameters == #
-n = 250
-replications = 100000
-distribution = uniform(loc=0, scale=(np.pi / 2))
-mu, s = distribution.mean(), distribution.std()
-
-g = np.sin
-g_prime = np.cos
-
-# == Generate obs of sqrt{n} (g(\bar X_n) - g(\mu)) == #
-data = distribution.rvs((replications, n)) 
-sample_means = data.mean(axis=1)  # Compute mean of each row
-error_obs = np.sqrt(n) * (g(sample_means) - g(mu))
-
-# == Plot == #
-asymptotic_sd = g_prime(mu) * s
-fig, ax = plt.subplots()
-xmin = -3 * g_prime(mu) * s
-xmax = -xmin
-ax.set_xlim(xmin, xmax)
-ax.hist(error_obs, bins=60, alpha=0.5, normed=True)
-xgrid = np.linspace(xmin, xmax, 200)
-lb = r"$N(0, g'(\mu)^2  \sigma^2)$"
-ax.plot(xgrid, norm.pdf(xgrid, scale=asymptotic_sd), 'k-', lw=2, label=lb)
-ax.legend()
-
-plt.show()
-
-
-
diff --git a/solutions/stand_alone_programs/solution_lln_ex2.py b/solutions/stand_alone_programs/solution_lln_ex2.py
deleted file mode 100644
index 308a662db..000000000
--- a/solutions/stand_alone_programs/solution_lln_ex2.py
+++ /dev/null
@@ -1,51 +0,0 @@
-"""
-Illustrates a consequence of the vector CLT.  The underlying random vector is
-X = (W, U + W), where W is Uniform(-1, 1), U is Uniform(-2, 2), and U and W
-are independent of each other.
-"""
-import numpy as np
-from scipy.stats import uniform, chi2
-from scipy.linalg import inv, sqrtm
-import matplotlib.pyplot as plt
-
-# == Set parameters == #
-n = 250
-replications = 50000
-dw = uniform(loc=-1, scale=2)  # Uniform(-1, 1)
-du = uniform(loc=-2, scale=4)  # Uniform(-2, 2)
-sw, su = dw.std(), du.std()
-vw, vu = sw**2, su**2
-Sigma = ((vw, vw), (vw, vw + vu))
-Sigma = np.array(Sigma)
-
-# == Compute Sigma^{-1/2} == #
-Q = inv(sqrtm(Sigma))  
-
-# == Generate observations of the normalized sample mean == #
-error_obs = np.empty((2, replications))
-for i in range(replications):
-    # == Generate one sequence of bivariate shocks == #
-    X = np.empty((2, n))
-    W = dw.rvs(n)
-    U = du.rvs(n)
-    # == Construct the n observations of the random vector == #
-    X[0, :] = W
-    X[1, :] = W + U
-    # == Construct the i-th observation of Y_n == #
-    error_obs[:, i] = np.sqrt(n) * X.mean(axis=1)
-
-# == Premultiply by Q and then take the squared norm == #
-temp = np.dot(Q, error_obs)
-chisq_obs = np.sum(temp**2, axis=0)
-
-# == Plot == #
-fig, ax = plt.subplots()
-xmax = 8
-ax.set_xlim(0, xmax)
-xgrid = np.linspace(0, xmax, 200)
-lb = "Chi-squared with 2 degrees of freedom"
-ax.plot(xgrid, chi2.pdf(xgrid, 2), 'k-', lw=2, label=lb)
-ax.legend()
-ax.hist(chisq_obs, bins=50, normed=True)
-
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_lqc_ex1.py b/solutions/stand_alone_programs/solution_lqc_ex1.py
deleted file mode 100644
index c93ceaeee..000000000
--- a/solutions/stand_alone_programs/solution_lqc_ex1.py
+++ /dev/null
@@ -1,80 +0,0 @@
-"""
-An LQ permanent income / life-cycle model with hump-shaped income
-
-    y_t = m1 * t + m2 * t^2 + sigma w_{t+1}
-
-where {w_t} is iid N(0, 1) and the coefficients m1 and m2 are chosen so that
-p(t) = m1 * t + m2 * t^2 has an inverted U shape with p(0) = 0, p(T/2) = mu
-and p(T) = 0.
-"""
-
-from __future__ import division
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lqcontrol import *
-
-# == Model parameters == #
-r       = 0.05
-beta    = 1 / (1 + r)
-T       = 50
-c_bar   = 1.5
-sigma   = 0.15
-mu      = 2
-q       = 1e4
-m1      = T * (mu / (T/2)**2)
-m2      = - (mu / (T/2)**2)
-
-# == Formulate as an LQ problem == #
-Q = 1
-R = np.zeros((4, 4)) 
-Rf = np.zeros((4, 4))
-Rf[0, 0] = q
-A = [[1 + r, -c_bar, m1, m2], 
-     [0,     1,      0,  0],
-     [0,     1,      1,  0],
-     [0,     1,      2,  1]]
-B = [[-1],
-     [0],
-     [0],
-     [0]]
-C = [[sigma],
-     [0],
-     [0],
-     [0]]
-
-# == Compute solutions and simulate == #
-lq = LQ(Q, R, A, B, C, beta=beta, T=T, Rf=Rf)
-x0 = (0, 1, 0, 0)
-xp, up, wp = lq.compute_sequence(x0)
-
-# == Convert results back to assets, consumption and income == #
-ap = xp[0, :]               # Assets
-c = up.flatten() + c_bar    # Consumption
-time = np.arange(1, T+1)
-income = wp[0, 1:] + m1 * time + m2 * time**2  # Income
-
-
-# == Plot results == #
-n_rows = 2
-fig, axes = plt.subplots(n_rows, 1, figsize=(12, 10))
-
-plt.subplots_adjust(hspace=0.5)
-for i in range(n_rows):
-    axes[i].grid()
-    axes[i].set_xlabel(r'Time')
-bbox = (0., 1.02, 1., .102)
-legend_args = {'bbox_to_anchor' : bbox, 'loc' : 3, 'mode' : 'expand'}
-p_args = {'lw' : 2, 'alpha' : 0.7}
-
-axes[0].plot(range(1, T+1), income, 'g-', label="non-financial income", **p_args)
-axes[0].plot(range(T), c, 'k-', label="consumption", **p_args)
-axes[1].plot(range(T+1), np.zeros(T+1), 'k-')
-axes[0].legend(ncol=2, **legend_args)
-
-axes[1].plot(range(T+1), ap.flatten(), 'b-', label="assets", **p_args)
-axes[1].plot(range(T), np.zeros(T), 'k-')
-axes[1].legend(ncol=1, **legend_args)
-
-plt.show()
-
-
diff --git a/solutions/stand_alone_programs/solution_lqc_ex2.py b/solutions/stand_alone_programs/solution_lqc_ex2.py
deleted file mode 100644
index 451dc0c1a..000000000
--- a/solutions/stand_alone_programs/solution_lqc_ex2.py
+++ /dev/null
@@ -1,107 +0,0 @@
-"""
-An permanent income / life-cycle model with polynomial growth in income
-over working life followed by a fixed retirement income.  The model is solved
-by combining two LQ programming problems as described in the lecture.
-"""
-
-from __future__ import division
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lqcontrol import *
-
-# == Model parameters == #
-r       = 0.05
-beta    = 1 / (1 + r)
-T       = 60
-K       = 40
-c_bar   = 4
-sigma   = 0.35
-mu      = 4
-q       = 1e4
-s       = 1
-m1      = 2 * mu / K
-m2      = - mu / K**2
-
-# == Formulate LQ problem 1 (retirement) == #
-Q = 1
-R = np.zeros((4, 4)) 
-Rf = np.zeros((4, 4))
-Rf[0, 0] = q
-A = [[1 + r, s - c_bar, 0, 0], 
-     [0,     1,      0,  0],
-     [0,     1,      1,  0],
-     [0,     1,      2,  1]]
-B = [[-1],
-     [0],
-     [0],
-     [0]]
-C = [[0],
-     [0],
-     [0],
-     [0]]
-
-# == Initialize LQ instance for retired agent == #
-lq_retired = LQ(Q, R, A, B, C, beta=beta, T=T-K, Rf=Rf)
-# == Iterate back to start of retirement, record final value function == #
-for i in range(T-K):
-    lq_retired.update_values()
-Rf2 = lq_retired.P
-
-# == Formulate LQ problem 2 (working life) == #
-R = np.zeros((4, 4)) 
-A = [[1 + r, -c_bar, m1, m2], 
-     [0,     1,      0,  0],
-     [0,     1,      1,  0],
-     [0,     1,      2,  1]]
-B = [[-1],
-     [0],
-     [0],
-     [0]]
-C = [[sigma],
-     [0],
-     [0],
-     [0]]
-
-# == Set up working life LQ instance with terminal Rf from lq_retired == #
-lq_working = LQ(Q, R, A, B, C, beta=beta, T=K, Rf=Rf2)
-
-# == Simulate working state / control paths == #
-x0 = (0, 1, 0, 0)
-xp_w, up_w, wp_w = lq_working.compute_sequence(x0)
-# == Simulate retirement paths (note the initial condition) == #
-xp_r, up_r, wp_r = lq_retired.compute_sequence(xp_w[:, K]) 
-
-# == Convert results back to assets, consumption and income == #
-xp = np.column_stack((xp_w, xp_r[:, 1:]))
-assets = xp[0, :]               # Assets
-
-up = np.column_stack((up_w, up_r))
-c = up.flatten() + c_bar    # Consumption
-
-time = np.arange(1, K+1)
-income_w = wp_w[0, 1:K+1] + m1 * time + m2 * time**2  # Income
-income_r = np.ones(T-K) * s
-income = np.concatenate((income_w, income_r))
-
-# == Plot results == #
-n_rows = 2
-fig, axes = plt.subplots(n_rows, 1, figsize=(12, 10))
-
-plt.subplots_adjust(hspace=0.5)
-for i in range(n_rows):
-    axes[i].grid()
-    axes[i].set_xlabel(r'Time')
-bbox = (0., 1.02, 1., .102)
-legend_args = {'bbox_to_anchor' : bbox, 'loc' : 3, 'mode' : 'expand'}
-p_args = {'lw' : 2, 'alpha' : 0.7}
-
-axes[0].plot(range(1, T+1), income, 'g-', label="non-financial income", **p_args)
-axes[0].plot(range(T), c, 'k-', label="consumption", **p_args)
-axes[1].plot(range(T+1), np.zeros(T+1), 'k-')
-axes[0].legend(ncol=2, **legend_args)
-
-axes[1].plot(range(T+1), assets, 'b-', label="assets", **p_args)
-axes[1].plot(range(T), np.zeros(T), 'k-')
-axes[1].legend(ncol=1, **legend_args)
-
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_lqc_ex3.py b/solutions/stand_alone_programs/solution_lqc_ex3.py
deleted file mode 100644
index fac9a39a6..000000000
--- a/solutions/stand_alone_programs/solution_lqc_ex3.py
+++ /dev/null
@@ -1,68 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: solution_lqc_ex3.py
-An infinite horizon profit maximization problem for a monopolist with
-adjustment costs.
-"""
-
-from __future__ import division
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lqcontrol import *
-
-# == Model parameters == #
-a0      = 5
-a1      = 0.5
-sigma   = 0.15
-rho     = 0.9
-gamma   = 1
-beta    = 0.95
-c       = 2
-T       = 120
-
-# == Useful constants == #
-m0 = (a0 - c) / (2 * a1)
-m1 = 1 / (2 * a1)
-
-# == Formulate LQ problem == #
-Q = gamma
-R = [[a1, -a1, 0],
-     [-a1, a1, 0],
-     [0,   0,  0]]
-A = [[rho, 0, m0 * (1 - rho)],
-     [0,   1, 0],
-     [0,   0, 1]]
-
-B = [[0],
-     [1],
-     [0]]
-C = [[m1 * sigma],
-     [0],
-     [0]]
-
-lq = LQ(Q, R, A, B, C=C, beta=beta)
-
-# == Simulate state / control paths == #
-x0 = (m0, 2, 1)
-xp, up, wp = lq.compute_sequence(x0, ts_length=150)
-q_bar = xp[0, :] 
-q     = xp[1, :]
-
-# == Plot simulation results == #
-fig, ax = plt.subplots(figsize=(10, 6.5))
-ax.set_xlabel('Time')
-
-# == Some fancy plotting stuff -- simplify if you prefer == #
-bbox = (0., 1.01, 1., .101)
-legend_args = {'bbox_to_anchor' : bbox, 'loc' : 3, 'mode' : 'expand'}
-p_args = {'lw' : 2, 'alpha' : 0.6}
-
-time = range(len(q))
-ax.set_xlim(0, max(time))
-ax.plot(time, q_bar, 'k-', lw=2, alpha=0.6, label=r'$\bar q_t$')
-ax.plot(time, q, 'b-', lw=2, alpha=0.6, label=r'$q_t$')
-ax.legend(ncol=2, **legend_args)
-s = r'dynamics with $\gamma = {}$'.format(gamma)
-ax.text(max(time) * 0.6, 1 * q_bar.max(), s, fontsize=14)
-
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_lqramsey_ex1.py b/solutions/stand_alone_programs/solution_lqramsey_ex1.py
deleted file mode 100644
index e627f5919..000000000
--- a/solutions/stand_alone_programs/solution_lqramsey_ex1.py
+++ /dev/null
@@ -1,33 +0,0 @@
-
-import numpy as np
-from numpy import array
-from quantecon.lqramsey import *
-
-# == Parameters == #
-beta = 1 / 1.05   
-rho, mg = .95, .35
-A = array([[0, 0, 0, rho, mg*(1-rho)],
-           [1, 0, 0, 0, 0],
-           [0, 1, 0, 0, 0],
-           [0, 0, 1, 0, 0],
-           [0, 0, 0, 0, 1]])
-C = np.zeros((5, 1))
-C[0, 0] = np.sqrt(1 - rho**2) * mg / 8
-Sg = array((1, 0, 0, 0, 0)).reshape(1, 5)        
-Sd = array((0, 0, 0, 0, 0)).reshape(1, 5)       
-Sb = array((0, 0, 0, 0, 2.135)).reshape(1, 5)  # Chosen st. (Sc + Sg) * x0 = 1
-Ss = array((0, 0, 0, 0, 0)).reshape(1, 5)
-
-economy = Economy(beta=beta, 
-        Sg=Sg, 
-        Sd=Sd, 
-        Sb=Sb, 
-        Ss=Ss, 
-        discrete=False, 
-        proc=(A, C))
-
-T = 50
-path = compute_paths(T, economy)
-gen_fig_1(path)
-
-
diff --git a/solutions/stand_alone_programs/solution_lss_ex1.py b/solutions/stand_alone_programs/solution_lss_ex1.py
deleted file mode 100644
index 5f85122e7..000000000
--- a/solutions/stand_alone_programs/solution_lss_ex1.py
+++ /dev/null
@@ -1,23 +0,0 @@
-
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lss import LSS
-
-phi_0, phi_1, phi_2 = 1.1, 0.8, -0.8
-
-A = [[1,     0,     0],
-     [phi_0, phi_1, phi_2],
-     [0,     1,     0]]
-C = np.zeros((3, 1))
-G = [0, 1, 0]
-
-ar = LSS(A, C, G, mu_0=np.ones(3))
-x, y = ar.simulate(ts_length=50)
-
-fig, ax = plt.subplots(figsize=(8, 4.6))
-y = y.flatten()
-ax.plot(y, 'b-', lw=2, alpha=0.7)
-ax.grid()
-ax.set_xlabel('time')
-ax.set_ylabel(r'$y_t$', fontsize=16)
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_lss_ex2.py b/solutions/stand_alone_programs/solution_lss_ex2.py
deleted file mode 100644
index 287eb37e0..000000000
--- a/solutions/stand_alone_programs/solution_lss_ex2.py
+++ /dev/null
@@ -1,29 +0,0 @@
-
-
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lss import LSS
-
-phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
-sigma = 0.2
-
-A = [[phi_1, phi_2, phi_3, phi_4],
-     [1,     0,     0,     0],
-     [0,     1,     0,     0],
-     [0,     0,     1,     0]]
-C = [[sigma], 
-     [0], 
-     [0], 
-     [0]]
-G = [1, 0, 0, 0]
-
-ar = LSS(A, C, G, mu_0=np.ones(4))
-x, y = ar.simulate(ts_length=200)
-
-fig, ax = plt.subplots(figsize=(8, 4.6))
-y = y.flatten()
-ax.plot(y, 'b-', lw=2, alpha=0.7)
-ax.grid()
-ax.set_xlabel('time')
-ax.set_ylabel(r'$y_t$', fontsize=16)
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_lss_ex3.py b/solutions/stand_alone_programs/solution_lss_ex3.py
deleted file mode 100644
index 2efdb8969..000000000
--- a/solutions/stand_alone_programs/solution_lss_ex3.py
+++ /dev/null
@@ -1,51 +0,0 @@
-
-from __future__ import division
-import numpy as np
-import matplotlib.pyplot as plt
-from scipy.stats import norm
-from quantecon.lss import LSS
-import random
-
-phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
-sigma = 0.1
-
-A = [[phi_1, phi_2, phi_3, phi_4],
-     [1,     0,     0,     0],
-     [0,     1,     0,     0],
-     [0,     0,     1,     0]]
-C = [[sigma], 
-     [0], 
-     [0], 
-     [0]]
-G = [1, 0, 0, 0]
-
-I = 20
-T = 50
-ar = LSS(A, C, G, mu_0=np.ones(4))
-ymin, ymax = -0.5, 1.15
-
-fig, ax = plt.subplots()
-
-ax.set_ylim(ymin, ymax)
-ax.set_xlabel(r'time', fontsize=16)
-ax.set_ylabel(r'$y_t$', fontsize=16)
-
-ensemble_mean = np.zeros(T)
-for i in range(I):
-    x, y = ar.simulate(ts_length=T)
-    y = y.flatten()
-    ax.plot(y, 'c-', lw=0.8, alpha=0.5)
-    ensemble_mean = ensemble_mean + y
-
-ensemble_mean = ensemble_mean / I
-ax.plot(ensemble_mean, color='b', lw=2, alpha=0.8, label=r'$\bar y_t$')
-
-m = ar.moment_sequence()
-population_means = []
-for t in range(T):
-    mu_x, mu_y, Sigma_x, Sigma_y = m.next()
-    population_means.append(mu_y)
-ax.plot(population_means, color='g', lw=2, alpha=0.8, label=r'$G\mu_t$')
-ax.legend(ncol=2)
-
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_lss_ex4.py b/solutions/stand_alone_programs/solution_lss_ex4.py
deleted file mode 100644
index bef00db28..000000000
--- a/solutions/stand_alone_programs/solution_lss_ex4.py
+++ /dev/null
@@ -1,49 +0,0 @@
-
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lss import LSS
-import random
-
-phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
-sigma = 0.1
-
-A = [[phi_1, phi_2, phi_3, phi_4],
-     [1,     0,     0,     0],
-     [0,     1,     0,     0],
-     [0,     0,     1,     0]]
-C = [[sigma], 
-     [0], 
-     [0], 
-     [0]]
-G = [1, 0, 0, 0]
-
-T0 = 10
-T1 = 50
-T2 = 75
-T4 = 100
-
-ar = LSS(A, C, G, mu_0=np.ones(4))
-ymin, ymax = -0.6, 0.6
-
-fig, ax = plt.subplots(figsize=(8, 5))
-
-ax.grid(alpha=0.4)
-ax.set_ylim(ymin, ymax)
-ax.set_ylabel(r'$y_t$', fontsize=16)
-ax.vlines((T0, T1, T2), -1.5, 1.5)
-
-ax.set_xticks((T0, T1, T2))
-ax.set_xticklabels((r"$T$", r"$T'$", r"$T''$"), fontsize=14)
-
-mu_x, mu_y, Sigma_x, Sigma_y = ar.stationary_distributions()
-ar.mu_0 = mu_x
-ar.Sigma_0 = Sigma_x
-
-for i in range(80):
-    rcolor = random.choice(('c', 'g', 'b'))
-    x, y = ar.simulate(ts_length=T4)
-    y = y.flatten()
-    ax.plot(y, color=rcolor, lw=0.8, alpha=0.5)
-    ax.plot((T0, T1, T2), (y[T0], y[T1], y[T2],), 'ko', alpha=0.5)
-
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_mass_ex1.py b/solutions/stand_alone_programs/solution_mass_ex1.py
deleted file mode 100644
index e281a8fce..000000000
--- a/solutions/stand_alone_programs/solution_mass_ex1.py
+++ /dev/null
@@ -1,32 +0,0 @@
-"""
-Filename: solution_mass_ex1.py
-Authors: David Evans, John Stachurski and Thomas J. Sargent
-LastModified: 12/02/2014
-"""
-
-import numpy as np
-from quantecon import asset_pricing 
-
-# == Define primitives == #
-n = 5
-P = 0.0125 * np.ones((n, n))
-P += np.diag(0.95 - 0.0125 * np.ones(5))
-s = np.array([1.05, 1.025, 1.0, 0.975, 0.95])
-gamma = 2.0
-beta = 0.94
-zeta = 1.0
-
-ap = asset_pricing.AssetPrices(beta, P, s, gamma)
-
-v = ap.tree_price()
-print "Lucas Tree Prices: ", v
-
-v_consol = ap.consol_price(zeta)
-print "Consol Bond Prices: ", v_consol
-
-P_tilde = P * s**(1-gamma)
-temp = beta * P_tilde.dot(v) - beta * P_tilde.dot(np.ones(n))
-print "Should be 0: ",  v - temp 
-
-p_s = 150.0
-w_bar, w_bars = ap.call_option(zeta, p_s, T = [10,20,30])
diff --git a/solutions/stand_alone_programs/solution_mass_ex2.py b/solutions/stand_alone_programs/solution_mass_ex2.py
deleted file mode 100644
index 6bc78174a..000000000
--- a/solutions/stand_alone_programs/solution_mass_ex2.py
+++ /dev/null
@@ -1,19 +0,0 @@
-from __future__ import division  # Omit for Python 3.x
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.lucastree import lucas_tree, compute_price
-
-fig, ax = plt.subplots()
-
-ax.set_xlabel(r'$y$', fontsize=16)
-ax.set_ylabel(r'price', fontsize=16)
-
-for beta in (.95, 0.98):
-    tree = lucas_tree(gamma=2, beta=beta, alpha=0.90, sigma=0.1)
-    grid, price_vals = compute_price(tree)
-    label = r'$\beta = {}$'.format(beta)
-    ax.plot(grid, price_vals, lw=2, alpha=0.7, label=label)
-
-ax.legend(loc='upper left')
-ax.set_xlim(min(grid), max(grid))
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_mc_ex1.py b/solutions/stand_alone_programs/solution_mc_ex1.py
deleted file mode 100644
index af081a622..000000000
--- a/solutions/stand_alone_programs/solution_mc_ex1.py
+++ /dev/null
@@ -1,33 +0,0 @@
-"""
-Compute the fraction of time that the worker spends unemployed,
-and compare it to the stationary probability.
-"""
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon import mc_tools
-
-alpha = beta = 0.1
-N = 10000
-p = beta / (alpha + beta)
-
-P = ((1 - alpha, alpha),   # Careful: P and p are distinct
-     (beta, 1 - beta))
-P = np.array(P)
-
-fig, ax = plt.subplots()
-ax.set_ylim(-0.25, 0.25)
-ax.grid()
-ax.hlines(0, 0, N, lw=2, alpha=0.6)  # Horizonal line at zero
-
-for x0, col in ((0, 'blue'), (1, 'green')):
-    # == Generate time series for worker that starts at x0 == #
-    X = mc_tools.sample_path(P, x0, N)
-    # == Compute fraction of time spent unemployed, for each n == #
-    X_bar = (X == 0).cumsum() / (1 + np.arange(N, dtype=float)) 
-    # == Plot == #
-    ax.fill_between(range(N), np.zeros(N), X_bar - p, color=col, alpha=0.1)
-    ax.plot(X_bar - p, color=col, label=r'$X_0 = \, {} $'.format(x0))
-    ax.plot(X_bar - p, 'k-', alpha=0.6)  # Overlay in black--make lines clearer
-
-ax.legend(loc='upper right')
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_mc_ex2.py b/solutions/stand_alone_programs/solution_mc_ex2.py
deleted file mode 100644
index 3fabbefe2..000000000
--- a/solutions/stand_alone_programs/solution_mc_ex2.py
+++ /dev/null
@@ -1,37 +0,0 @@
-"""
-Return list of pages, ordered by rank
-"""
-from __future__ import print_function, division  # Omit if using Python 3.x
-import numpy as np
-from quantecon import mc_tools
-from operator import itemgetter
-import re
-
-infile = 'web_graph_data.txt'
-alphabet = 'abcdefghijklmnopqrstuvwxyz'
-
-n = 14 # Total number of web pages (nodes)
-
-# == Create a matrix Q indicating existence of links == #
-#  * Q[i, j] = 1 if there is a link from i to j
-#  * Q[i, j] = 0 otherwise
-Q = np.zeros((n, n), dtype=int)
-f = open(infile, 'r')
-edges = f.readlines()
-f.close()
-for edge in edges:
-    from_node, to_node = re.findall('\w', edge)
-    i, j = alphabet.index(from_node), alphabet.index(to_node)
-    Q[i, j] = 1
-# == Create the corresponding Markov matrix P == #
-P = np.empty((n, n))
-for i in range(n):
-    P[i,:] = Q[i,:] / Q[i,:].sum()
-# == Compute the stationary distribution r == #
-r = mc_tools.compute_stationary(P)
-ranked_pages = {alphabet[i] : r[i] for i in range(n)}
-# == Print solution, sorted from highest to lowest rank == #
-print('Rankings\n ***')
-for name, rank in sorted(ranked_pages.iteritems(), key=itemgetter(1), reverse=1):
-    print('{0}: {1:.4}'.format(name, rank))
-
diff --git a/solutions/stand_alone_programs/solution_odu_ex1.py b/solutions/stand_alone_programs/solution_odu_ex1.py
deleted file mode 100644
index 4e5e4f49e..000000000
--- a/solutions/stand_alone_programs/solution_odu_ex1.py
+++ /dev/null
@@ -1,50 +0,0 @@
-"""
-Solves the "Offer Distribution Unknown" model by iterating on a guess of the
-reservation wage function.
-"""
-from scipy import interp
-import numpy as np
-from numpy import maximum as npmax
-import matplotlib.pyplot as plt
-from quantecon.odu_vfi import searchProblem
-from scipy.integrate import fixed_quad
-from quantecon.compute_fp import compute_fixed_point
-
-
-def res_wage_operator(sp, phi):
-    """
-    Updates the reservation wage function guess phi via the operator Q.
-    Returns the updated function Q phi, represented as the array new_phi.
-    
-        * sp is an instance of searchProblem, defined in odu_vfi
-        * phi is a NumPy array with len(phi) = len(sp.pi_grid)
-
-    """
-    beta, c, f, g, q = sp.beta, sp.c, sp.f, sp.g, sp.q    # Simplify names
-    phi_f = lambda p: interp(p, sp.pi_grid, phi)  # Turn phi into a function
-    new_phi = np.empty(len(phi))
-    for i, pi in enumerate(sp.pi_grid):
-        def integrand(x):
-            "Integral expression on right-hand side of operator"
-            return npmax(x, phi_f(q(x, pi))) * (pi * f(x) + (1 - pi) * g(x))
-        integral, error = fixed_quad(integrand, 0, sp.w_max)
-        new_phi[i] = (1 - beta) * c + beta * integral
-    return new_phi
-
-
-if __name__ == '__main__':  # If module is run rather than imported
-
-    sp = searchProblem(pi_grid_size=50)
-    phi_init = np.ones(len(sp.pi_grid)) 
-    w_bar = compute_fixed_point(res_wage_operator, sp, phi_init)
-
-    fig, ax = plt.subplots()
-    ax.plot(sp.pi_grid, w_bar, linewidth=2, color='black')
-    ax.set_ylim(0, 2)
-    ax.grid(axis='x', linewidth=0.25, linestyle='--', color='0.25')
-    ax.grid(axis='y', linewidth=0.25, linestyle='--', color='0.25')
-    ax.fill_between(sp.pi_grid, 0, w_bar, color='blue', alpha=0.15)
-    ax.fill_between(sp.pi_grid, w_bar, 2, color='green', alpha=0.15)
-    ax.text(0.42, 1.2, 'reject')
-    ax.text(0.7, 1.8, 'accept')
-    plt.show()
diff --git a/solutions/stand_alone_programs/solution_odu_ex2.py b/solutions/stand_alone_programs/solution_odu_ex2.py
deleted file mode 100644
index cbb4d9305..000000000
--- a/solutions/stand_alone_programs/solution_odu_ex2.py
+++ /dev/null
@@ -1,65 +0,0 @@
-from scipy import interp
-import numpy as np
-import matplotlib.pyplot as plt
-from quantecon.odu_vfi import searchProblem
-from quantecon.compute_fp import compute_fixed_point
-from solution_odu_ex1 import res_wage_operator
-
-
-# Set up model and compute the function w_bar
-sp = searchProblem(pi_grid_size=50, F_a=1, F_b=1)
-pi_grid, f, g, F, G = sp.pi_grid, sp.f, sp.g, sp.F, sp.G
-phi_init = np.ones(len(sp.pi_grid)) 
-w_bar_vals = compute_fixed_point(res_wage_operator, sp, phi_init)
-w_bar = lambda x: interp(x, pi_grid, w_bar_vals)
-
-
-class Agent:
-    """
-    Holds the employment state and beliefs of an individual agent.
-    """
-
-    def __init__(self, pi=1e-3):
-        self.pi = pi
-        self.employed = 1
-
-    def update(self, H):
-        "Update self by drawing wage offer from distribution H."
-        if self.employed == 0:
-            w = H.rvs()
-            if w >= w_bar(self.pi):
-                self.employed = 1
-            else:
-                self.pi = 1.0 / (1 + ((1 - self.pi) * g(w)) / (self.pi * f(w)))
-
-
-num_agents = 5000
-separation_rate = 0.025  # Fraction of jobs that end in each period 
-separation_num = int(num_agents * separation_rate)
-agent_indices = range(num_agents)
-agents = [Agent() for i in range(num_agents)]
-sim_length = 600
-H = G  # Start with distribution G
-change_date = 200  # Change to F after this many periods
-
-unempl_rate = []
-for i in range(sim_length):
-    print "date = ", i
-    if i == change_date:
-        H = F
-    # Randomly select separation_num agents and set employment status to 0
-    np.random.shuffle(agent_indices)
-    separation_list = agent_indices[:separation_num]
-    for agent_index in separation_list:
-        agents[agent_index].employed = 0
-    # Update agents
-    for agent in agents:
-        agent.update(H)
-    employed = [agent.employed for agent in agents]
-    unempl_rate.append(1 - np.mean(employed))
-
-fig, ax = plt.subplots()
-ax.plot(unempl_rate, lw=2, alpha=0.8, label='unemployment rate')
-ax.axvline(change_date, color="red")
-ax.legend()
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_og_ex1.py b/solutions/stand_alone_programs/solution_og_ex1.py
deleted file mode 100644
index 7d61d47ec..000000000
--- a/solutions/stand_alone_programs/solution_og_ex1.py
+++ /dev/null
@@ -1,26 +0,0 @@
-import matplotlib.pyplot as plt
-from quantecon.optgrowth import growthModel, bellman_operator, compute_greedy
-from quantecon.compute_fp import compute_fixed_point
-
-alpha, beta = 0.65, 0.95
-gm = growthModel() 
-true_sigma = (1 - alpha * beta) * gm.grid**alpha
-w = 5 * gm.u(gm.grid) - 25  # Initial condition
-
-fig, ax = plt.subplots(3, 1, figsize=(8, 10))
-
-for i, n in enumerate((2, 4, 6)):
-    ax[i].set_ylim(0, 1)
-    ax[i].set_xlim(0, 2)
-    ax[i].set_yticks((0, 1))
-    ax[i].set_xticks((0, 2))
-
-    v_star = compute_fixed_point(bellman_operator, gm, w, max_iter=n)
-    sigma = compute_greedy(gm, v_star)
-
-    ax[i].plot(gm.grid, sigma, 'b-', lw=2, alpha=0.8, label='approximate optimal policy')
-    ax[i].plot(gm.grid, true_sigma, 'k-', lw=2, alpha=0.8, label='true optimal policy')
-    ax[i].legend(loc='upper left')
-    ax[i].set_title('{} value function iterations'.format(n))
-
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_og_ex2.py b/solutions/stand_alone_programs/solution_og_ex2.py
deleted file mode 100644
index 77f4a881a..000000000
--- a/solutions/stand_alone_programs/solution_og_ex2.py
+++ /dev/null
@@ -1,34 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-from scipy import interp
-from quantecon.optgrowth import growthModel, bellman_operator, compute_greedy
-from quantecon.compute_fp import compute_fixed_point
-
-gm = growthModel() 
-w = 5 * gm.u(gm.grid) - 25  # To be used as an initial condition
-discount_factors = (0.9, 0.94, 0.98)
-series_length = 25
-
-fig, ax = plt.subplots()
-ax.set_xlabel("time")
-ax.set_ylabel("capital")
-
-for beta in discount_factors:
-
-    # Compute the optimal policy given the discount factor
-    gm.beta = beta
-    v_star = compute_fixed_point(bellman_operator, gm, w, max_iter=20)
-    sigma = compute_greedy(gm, v_star)
-
-    # Compute the corresponding time series for capital
-    k = np.empty(series_length)
-    k[0] = 0.1
-    sigma_function = lambda x: interp(x, gm.grid, sigma)
-    for t in range(1, series_length):
-        k[t] = gm.f(k[t-1]) - sigma_function(k[t-1])
-    ax.plot(k, 'o-', lw=2, alpha=0.75, label=r'$\beta = {}$'.format(beta))
-
-ax.legend(loc='lower right')
-plt.show()
-
-
diff --git a/solutions/stand_alone_programs/solution_oop_ex1.py b/solutions/stand_alone_programs/solution_oop_ex1.py
deleted file mode 100644
index 6239efe61..000000000
--- a/solutions/stand_alone_programs/solution_oop_ex1.py
+++ /dev/null
@@ -1,12 +0,0 @@
-class ecdf:
-
-    def __init__(self, observations):
-        self.observations = observations
-
-    def __call__(self, x):
-        counter = 0.0
-        for obs in self.observations:
-            if obs <= x:
-                counter += 1
-        return counter / len(self.observations)
-
diff --git a/solutions/stand_alone_programs/solution_oop_ex2.py b/solutions/stand_alone_programs/solution_oop_ex2.py
deleted file mode 100644
index d181988d1..000000000
--- a/solutions/stand_alone_programs/solution_oop_ex2.py
+++ /dev/null
@@ -1,29 +0,0 @@
-class Polynomial:
-
-    def __init__(self, coefficients):
-        """
-        Creates an instance of the Polynomial class representing 
-
-            p(x) = a_0 x^0 + ... + a_N x^N, 
-            
-        where a_i = coefficients[i].
-        """
-        self.coefficients = coefficients
-
-    def __call__(self, x):
-        "Evaluate the polynomial at x."
-        y = 0
-        for i, a in enumerate(self.coefficients):
-            y += a * x**i  
-        return y
-
-    def differentiate(self):
-        "Reset self.coefficients to those of p' instead of p."
-        new_coefficients = []
-        for i, a in enumerate(self.coefficients):
-            new_coefficients.append(i * a)
-        # Remove the first element, which is zero
-        del new_coefficients[0]  
-        # And reset coefficients data to new values
-        self.coefficients = new_coefficients
-
diff --git a/solutions/stand_alone_programs/solution_pbe_ex1.py b/solutions/stand_alone_programs/solution_pbe_ex1.py
deleted file mode 100644
index 7fc421fbf..000000000
--- a/solutions/stand_alone_programs/solution_pbe_ex1.py
+++ /dev/null
@@ -1,5 +0,0 @@
-def factorial(n):
-    k = 1
-    for i in range(n):
-        k = k * (i + 1)
-    return k
diff --git a/solutions/stand_alone_programs/solution_pbe_ex2.py b/solutions/stand_alone_programs/solution_pbe_ex2.py
deleted file mode 100644
index e6e246298..000000000
--- a/solutions/stand_alone_programs/solution_pbe_ex2.py
+++ /dev/null
@@ -1,10 +0,0 @@
-from random import uniform
-
-def binomial_rv(n, p):
-    count = 0
-    for i in range(n):
-        U = uniform(0, 1)
-        if U < p:
-            count = count + 1    # Or count += 1
-    print count
-
diff --git a/solutions/stand_alone_programs/solution_pbe_ex3.py b/solutions/stand_alone_programs/solution_pbe_ex3.py
deleted file mode 100644
index df0a45daa..000000000
--- a/solutions/stand_alone_programs/solution_pbe_ex3.py
+++ /dev/null
@@ -1,17 +0,0 @@
-from __future__ import division  # Omit if using Python 3.x
-from random import uniform
-from math import sqrt
-
-n = 100000
-
-count = 0
-for i in range(n):
-    u, v = uniform(0, 1), uniform(0, 1)
-    d = sqrt((u - 0.5)**2 + (v - 0.5)**2)
-    if d < 0.5:
-        count += 1
-
-area_estimate = count / n
-
-print area_estimate * 4  # dividing by radius**2
-
diff --git a/solutions/stand_alone_programs/solution_pbe_ex4.py b/solutions/stand_alone_programs/solution_pbe_ex4.py
deleted file mode 100644
index 87597a7ed..000000000
--- a/solutions/stand_alone_programs/solution_pbe_ex4.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from random import uniform
-
-payoff = 0
-count = 0
-
-for i in range(10):
-    U = uniform(0, 1)
-    count = count + 1 if U < 0.5 else 0
-    if count == 3:
-        payoff = 1
-
-print payoff
diff --git a/solutions/stand_alone_programs/solution_pbe_ex5.py b/solutions/stand_alone_programs/solution_pbe_ex5.py
deleted file mode 100644
index 0b333aadf..000000000
--- a/solutions/stand_alone_programs/solution_pbe_ex5.py
+++ /dev/null
@@ -1,14 +0,0 @@
-from pylab import plot, show
-from random import normalvariate
-
-alpha = 0.9
-ts_length = 200
-current_x = 0
-
-x_values = []
-for i in range(ts_length):
-    x_values.append(current_x)
-    current_x = alpha * current_x + normalvariate(0, 1)
-plot(x_values, 'b-')
-show()
-
diff --git a/solutions/stand_alone_programs/solution_pbe_ex6.py b/solutions/stand_alone_programs/solution_pbe_ex6.py
deleted file mode 100644
index dc4e5d435..000000000
--- a/solutions/stand_alone_programs/solution_pbe_ex6.py
+++ /dev/null
@@ -1,17 +0,0 @@
-from pylab import plot, show, legend
-from random import normalvariate
-
-alphas = [0.0, 0.8, 0.98]
-ts_length = 200
-
-for alpha in alphas:
-    x_values = []
-    current_x = 0
-    for i in range(ts_length):
-        x_values.append(current_x)
-        current_x = alpha * current_x + normalvariate(0, 1)
-    plot(x_values, label='alpha = ' + str(alpha))
-legend()
-show()
-
-
diff --git a/solutions/stand_alone_programs/solution_pd_ex1.py b/solutions/stand_alone_programs/solution_pd_ex1.py
deleted file mode 100644
index 607137098..000000000
--- a/solutions/stand_alone_programs/solution_pd_ex1.py
+++ /dev/null
@@ -1,39 +0,0 @@
-import numpy as np
-import pandas as pd
-import datetime as dt
-import pandas.io.data as web
-import matplotlib.pyplot as plt
-
-ticker_list = {'INTC': 'Intel',
-               'MSFT': 'Microsoft',
-               'IBM': 'IBM',
-               'BHP': 'BHP',
-               'RSH': 'RadioShack',
-               'TM': 'Toyota',
-               'AAPL': 'Apple',
-               'AMZN': 'Amazon',
-               'BA': 'Boeing',
-               'QCOM': 'Qualcomm',
-               'KO': 'Coca-Cola',
-               'GOOG': 'Google',
-               'SNE': 'Sony',
-               'PTR': 'PetroChina'}
-
-start = dt.datetime(2013, 1, 1)
-end = dt.datetime.today()
-
-price_change = {}
-
-for ticker in ticker_list:
-    prices = web.DataReader(ticker, 'yahoo', start, end)
-    closing_prices = prices['Close']
-    change = 100 * (closing_prices[-1] - closing_prices[0]) / closing_prices[0]
-    name = ticker_list[ticker]
-    price_change[name] = change
-
-pc = pd.Series(price_change)
-pc.sort()
-fig, ax = plt.subplots()
-pc.plot(kind='bar', ax=ax)
-plt.show()
-
diff --git a/solutions/stand_alone_programs/solution_ree_ex1.py b/solutions/stand_alone_programs/solution_ree_ex1.py
deleted file mode 100644
index e160e17bb..000000000
--- a/solutions/stand_alone_programs/solution_ree_ex1.py
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: solution_ree_ex1.py
-Authors: Chase Coleman, Spencer Lyon, Thomas Sargent, John Stachurski
-Solves an exercise from the rational expectations module
-"""
-from __future__ import print_function
-import numpy as np
-from quantecon.lqcontrol import LQ
-
-
-# == Model parameters == #
-
-a0      = 100
-a1      = 0.05
-beta    = 0.95
-gamma   = 10.0
-
-# == Beliefs == #
-
-kappa0  = 95.5
-kappa1  = 0.95
-
-# == Formulate the LQ problem == #
-
-A = np.array([[1, 0, 0], [0, kappa1, kappa0], [0, 0, 1]])
-B = np.array([1, 0, 0])
-B.shape = 3, 1
-R = np.array([[0, -a1/2, a0/2], [-a1/2, 0, 0], [a0/2, 0, 0]])
-Q = -0.5 * gamma
-
-# == Solve for the optimal policy == #
-
-lq = LQ(Q, R, A, B, beta=beta)
-P, F, d = lq.stationary_values()
-F = F.flatten()
-out1 = "F = [{0:.3f}, {1:.3f}, {2:.3f}]".format(F[0], F[1], F[2])
-h0, h1, h2 = -F[2], 1 - F[0], -F[1]
-out2 = "(h0, h1, h2) = ({0:.3f}, {1:.3f}, {2:.3f})".format(h0, h1, h2)
-
-print(out1)
-print(out2)
diff --git a/solutions/stand_alone_programs/solution_ree_ex2.py b/solutions/stand_alone_programs/solution_ree_ex2.py
deleted file mode 100644
index c7299c5c0..000000000
--- a/solutions/stand_alone_programs/solution_ree_ex2.py
+++ /dev/null
@@ -1,35 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: solution_ree_ex2.py
-Authors: Chase Coleman, Spencer Lyon, Thomas Sargent, John Stachurski
-Solves an exercise from the rational expectations module
-"""
-from __future__ import print_function
-import numpy as np
-from quantecon.lqcontrol import LQ
-from quantecon.solution_ree_ex1 import beta, R, Q, B
-
-candidates = (
-          (94.0886298678, 0.923409232937),
-          (93.2119845412, 0.984323478873),
-          (95.0818452486, 0.952459076301)
-             )
-
-for kappa0, kappa1 in candidates:
-
-    # == Form the associated law of motion == #
-    A = np.array([[1, 0, 0], [0, kappa1, kappa0], [0, 0, 1]])
-
-    # == Solve the LQ problem for the firm == #
-    lq = LQ(Q, R, A, B, beta=beta)
-    P, F, d = lq.stationary_values()
-    F = F.flatten()
-    h0, h1, h2 = -F[2], 1 - F[0], -F[1]
-
-    # == Test the equilibrium condition == #
-    if np.allclose((kappa0, kappa1), (h0, h1 + h2)):
-        print('Equilibrium pair =', kappa0, kappa1)
-        print('(h0, h1, h2) = ', h0, h1, h2)
-        break
-
-
diff --git a/solutions/stand_alone_programs/solution_ree_ex3.py b/solutions/stand_alone_programs/solution_ree_ex3.py
deleted file mode 100644
index 27e2fdf87..000000000
--- a/solutions/stand_alone_programs/solution_ree_ex3.py
+++ /dev/null
@@ -1,29 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: solution_ree_ex3.py
-Authors: Chase Coleman, Spencer Lyon, Thomas Sargent, John Stachurski
-Solves an exercise from the rational expectations module
-"""
-
-from __future__ import print_function
-import numpy as np
-from quantecon.lqcontrol import LQ
-from solution_ree_ex1 import a0, a1, beta, gamma
-
-# == Formulate the planner's LQ problem == #
-
-A = np.array([[1, 0], [0, 1]])
-B = np.array([[1], [0]])
-R = -np.array([[a1 / 2, -a0 / 2], [-a0 / 2, 0]])
-Q = - gamma / 2
-
-# == Solve for the optimal policy == #
-
-lq = LQ(Q, R, A, B, beta=beta)
-P, F, d = lq.stationary_values()
-
-# == Print the results == #
-
-F = F.flatten()
-kappa0, kappa1 = -F[1], 1 - F[0]
-print(kappa0, kappa1)
diff --git a/solutions/stand_alone_programs/solution_ree_ex4.py b/solutions/stand_alone_programs/solution_ree_ex4.py
deleted file mode 100644
index 4f1731a01..000000000
--- a/solutions/stand_alone_programs/solution_ree_ex4.py
+++ /dev/null
@@ -1,23 +0,0 @@
-"""
-Origin: QE by John Stachurski and Thomas J. Sargent
-Filename: solution_ree_ex4.py
-Authors: Chase Coleman, Spencer Lyon, Thomas Sargent, John Stachurski
-Solves an exercise from the rational expectations module
-"""
-
-from __future__ import print_function
-import numpy as np
-from quantecon.lqcontrol import LQ
-from solution_ree_ex1 import a0, a1, beta, gamma
-
-A = np.array([[1, 0], [0, 1]])
-B = np.array([[1], [0]])
-R = - np.array([[a1, -a0 / 2], [-a0 / 2, 0]])
-Q = - gamma / 2
-
-lq = LQ(Q, R, A, B, beta=beta)
-P, F, d = lq.stationary_values()
-
-F = F.flatten()
-m0, m1 = -F[1], 1 - F[0]
-print(m0, m1)
diff --git a/solutions/stand_alone_programs/solution_shortpath.py b/solutions/stand_alone_programs/solution_shortpath.py
deleted file mode 100644
index 63c21779d..000000000
--- a/solutions/stand_alone_programs/solution_shortpath.py
+++ /dev/null
@@ -1,79 +0,0 @@
-"""
-Source: QE by John Stachurski and Thomas J. Sargent
-Filename: solution_shortpath.py
-Authors: John Stachurksi and Thomas J. Sargent
-LastModified: 11/08/2013
-"""
-
-def read_graph():
-    """ Read in the graph from the data file.  The graph is stored
-    as a dictionary, where the keys are the nodes, and the values
-    are a list of pairs (d, c), where d is a node and c is a number.
-    If (d, c) is in the list for node n, then d can be reached from
-    n at cost c.
-    """
-    graph = {}
-    infile = open('graph.txt')
-    for line in infile:
-        elements = line.split(',')
-        node = elements.pop(0).strip()
-        graph[node] = []
-        if node != 'node99':
-            for element in elements:
-                destination, cost = element.split()
-                graph[node].append((destination.strip(), float(cost)))
-    infile.close()
-    return graph
-
-def update_J(J, graph):
-    "The Bellman operator."
-    next_J = {}
-    for node in graph:
-        if node == 'node99':
-            next_J[node] = 0
-        else:
-            next_J[node] = min(cost + J[dest] for dest, cost in graph[node])
-    return next_J
-
-def print_best_path(J, graph):
-    """ Given a cost-to-go function, computes the best path.  At each node n, 
-    the function prints the current location, looks at all nodes that can be 
-    reached from n, and moves to the node m which minimizes c + J[m], where c 
-    is the cost of moving to m.
-    """
-    sum_costs = 0
-    current_location = 'node0'
-    while current_location != 'node99':
-        print current_location
-        running_min = 1e100  # Any big number
-        for destination, cost in graph[current_location]:
-            cost_of_path = cost + J[destination]
-            if cost_of_path < running_min:
-                running_min = cost_of_path
-                minimizer_cost = cost
-                minimizer_dest = destination
-        current_location = minimizer_dest
-        sum_costs += minimizer_cost
-
-    print 'node99'
-    print
-    print 'Cost: ', sum_costs
-
-
-## Main loop
-
-graph = read_graph()
-M = 1e10
-J = {}
-for node in graph:
-    J[node] = M
-J['node99'] = 0
-
-while 1:
-    next_J = update_J(J, graph)
-    if next_J == J:
-        break
-    else:
-        J = next_J
-print_best_path(J, graph)
-
diff --git a/solutions/stand_alone_programs/solution_statd_ex1.py b/solutions/stand_alone_programs/solution_statd_ex1.py
deleted file mode 100644
index d6c9e7664..000000000
--- a/solutions/stand_alone_programs/solution_statd_ex1.py
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-Look ahead estimation of a TAR stationary density, where the TAR model is
-
-    X' = theta |X| + sqrt(1 - theta^2) xi
-
-and xi is standard normal.  Try running at n = 10, 100, 1000, 10000 to get an
-idea of the speed of convergence.
-"""
-import numpy as np
-from scipy.stats import norm, gaussian_kde
-import matplotlib.pyplot as plt
-from quantecon.lae import lae
-
-phi = norm()
-n = 500
-theta = 0.8
-# == Frequently used constants == #
-d = np.sqrt(1 - theta**2) 
-delta = theta / d
-
-def psi_star(y):
-    "True stationary density of the TAR Model"
-    return 2 * norm.pdf(y) * norm.cdf(delta * y) 
-
-def p(x, y):
-        "Stochastic kernel for the TAR model."
-        return phi.pdf((y - theta * np.abs(x)) / d) / d
-
-Z = phi.rvs(n)
-X = np.empty(n)
-for t in range(n-1):
-    X[t+1] = theta * np.abs(X[t]) + d * Z[t]
-psi_est = lae(p, X)
-k_est = gaussian_kde(X)
-
-fig, ax = plt.subplots()
-ys = np.linspace(-3, 3, 200)
-ax.plot(ys, psi_star(ys), 'b-', lw=2, alpha=0.6, label='true')
-ax.plot(ys, psi_est(ys), 'g-', lw=2, alpha=0.6, label='look ahead estimate')
-ax.plot(ys, k_est(ys), 'k-', lw=2, alpha=0.6, label='kernel based estimate')
-ax.legend(loc='upper left')
-plt.show()
diff --git a/solutions/stand_alone_programs/solution_statd_ex2.py b/solutions/stand_alone_programs/solution_statd_ex2.py
deleted file mode 100644
index e2dbde26f..000000000
--- a/solutions/stand_alone_programs/solution_statd_ex2.py
+++ /dev/null
@@ -1,51 +0,0 @@
-import numpy as np
-from scipy.stats import lognorm, beta
-import matplotlib.pyplot as plt
-from quantecon.lae import lae
-
-# == Define parameters == #
-s = 0.2
-delta = 0.1
-a_sigma = 0.4       # A = exp(B) where B ~ N(0, a_sigma)
-alpha = 0.4         # f(k) = k^{\alpha}
-
-phi = lognorm(a_sigma) 
-
-def p(x, y):
-    "Stochastic kernel, vectorized in x.  Both x and y must be positive."
-    d = s * x**alpha
-    return phi.pdf((y - (1 - delta) * x) / d) / d
-
-n = 1000     # Number of observations at each date t
-T = 40       # Compute density of k_t at 1,...,T
-
-fig, axes = plt.subplots(2, 2)
-axes = axes.flatten()
-xmax = 6.5
-
-for i in range(4):
-    ax = axes[i] 
-    ax.set_xlim(0, xmax)
-    psi_0 = beta(5, 5, scale=0.5, loc=i*2)  # Initial distribution
-
-    # == Generate matrix s.t. t-th column is n observations of k_t == #
-    k = np.empty((n, T))
-    A = phi.rvs((n, T))
-    k[:, 0] = psi_0.rvs(n)
-    for t in range(T-1):
-        k[:, t+1] = s * A[:,t] * k[:, t]**alpha + (1 - delta) * k[:, t]
-
-    # == Generate T instances of lae using this data, one for each t == #
-    laes = [lae(p, k[:, t]) for t in range(T)]
-
-    ygrid = np.linspace(0.01, xmax, 150)
-    greys = [str(g) for g in np.linspace(0.0, 0.8, T)]
-    greys.reverse()
-    for psi, g in zip(laes, greys):
-        ax.plot(ygrid, psi(ygrid), color=g, lw=2, alpha=0.6)
-    #ax.set_xlabel('capital')
-    #title = r'Density of $k_1$ (lighter) to $k_T$ (darker) for $T={}$'
-    #ax.set_title(title.format(T))
-
-plt.show()
-