Manual expectations

lectures/dynamic_programming/jv.md

@@ -44,7 +44,7 @@ kernelspec:
 tags: [hide-output]
 ---
 using LinearAlgebra, Statistics
-using Distributions, Interpolations, Expectations
+using Distributions, Interpolations
 using LaTeXStrings, Plots, NLsolve, Random
 
 ```
@@ -194,8 +194,9 @@ function JvWorker(; A = 1.4,
     pi_func = sqrt
     F = Beta(2, 2)
 
-    # expectation operator
-    E = expectation(F)
+    # probability vector and support for manual expectation
+    w = support(F)
+    p = pdf.(F, w)
 
     # Set up grid over the state space for DP
     # Max of grid is the max of a large quantile value for F and the
@@ -207,13 +208,13 @@ function JvWorker(; A = 1.4,
     x_grid = range(epsilon, grid_max, length = grid_size)
 
     return (; A, alpha, beta, x_grid, G,
-            pi_func, F, E, epsilon)
+            pi_func, F, w, p, epsilon)
 end
 
-function T!(jv, V, new_V::AbstractVector)
+function T!(jv, V, new_V)
 
     # simplify notation
-    (; G, pi_func, F, beta, E, epsilon) = jv
+    (; G, pi_func, beta, w, p, epsilon) = jv
 
     # prepare interpoland of value function
     Vf = LinearInterpolation(jv.x_grid, V, extrapolation_bc = Line())
@@ -225,19 +226,18 @@ function T!(jv, V, new_V::AbstractVector)
     max_phi = 1.0
     search_grid = range(epsilon, 1.0, length = 15)
 
-    for (i, x) in enumerate(jv.x_grid)
-        function w(z)
-            s, phi = z
-            h(u) = Vf(max(G(x, phi), u))
-            integral = E(h)
-            q = pi_func(s) * integral + (1.0 - pi_func(s)) * Vf(G(x, phi))
-
-            return -x * (1.0 - phi - s) - beta * q
-        end
+    # objective function
+    function w_x(x, s, phi)
+        h(u) = Vf(max(G(x, phi), u))
+        integral = dot(h.(w), p)
+        q = pi_func(s) * integral + (1.0 - pi_func(s)) * Vf(G(x, phi))
+        return -x * (1.0 - phi - s) - beta * q
+    end
 
+    for (i, x) in enumerate(jv.x_grid)
         for s in search_grid
             for phi in search_grid
-                cur_val = ifelse(s + phi <= 1.0, -w((s, phi)), -1.0)
+                cur_val = ifelse(s + phi <= 1.0, -w_x(x, s, phi), -1.0)
                 if cur_val > max_val
                     max_val, max_s, max_phi = cur_val, s, phi
                 end
@@ -248,10 +248,10 @@ function T!(jv, V, new_V::AbstractVector)
     end
 end
 
-function T!(jv, V, out::Tuple{AbstractVector, AbstractVector})
+function T!(jv, V, out::Tuple)
 
     # simplify notation
-    (; G, pi_func, F, beta, E, epsilon) = jv
+    (; G, pi_func, beta, w, p, epsilon) = jv
 
     # prepare interpoland of value function
     Vf = LinearInterpolation(jv.x_grid, V, extrapolation_bc = Line())
@@ -266,19 +266,18 @@ function T!(jv, V, out::Tuple{AbstractVector, AbstractVector})
     max_phi = 1.0
     search_grid = range(epsilon, 1.0, length = 15)
 
-    for (i, x) in enumerate(jv.x_grid)
-        function w(z)
-            s, phi = z
-            h(u) = Vf(max(G(x, phi), u))
-            integral = E(h)
-            q = pi_func(s) * integral + (1.0 - pi_func(s)) * Vf(G(x, phi))
-
-            return -x * (1.0 - phi - s) - beta * q
-        end
+    # objective function
+    function w_x(x, s, phi)
+        h(u) = Vf(max(G(x, phi), u))
+        integral = dot(h.(w), p)
+        q = pi_func(s) * integral + (1.0 - pi_func(s)) * Vf(G(x, phi))
+        return -x * (1.0 - phi - s) - beta * q
+    end
 
+    for (i, x) in enumerate(jv.x_grid)
         for s in search_grid
             for phi in search_grid
-                cur_val = ifelse(s + phi <= 1.0, -w((s, phi)), -1.0)
+                cur_val = ifelse(s + phi <= 1.0, -w_x(x, s, phi), -1.0)
                 if cur_val > max_val
                     max_val, max_s, max_phi = cur_val, s, phi
                 end

lectures/dynamic_programming/mccall_model.md

@@ -291,7 +291,7 @@ generates a sequence that converges to the fixed point.
 tags: [hide-output]
 ---
 using LinearAlgebra, Statistics
-using Distributions, Expectations, LaTeXStrings
+using Distributions, LaTeXStrings
 using NLsolve, Roots, Random, Plots
 ```
 
@@ -312,27 +312,26 @@ Here's the distribution of wage offers we'll work with
 n = 50
 dist = BetaBinomial(n, 200, 100) # probability distribution
 @show support(dist)
+p = pdf.(dist, support(dist))
 w = range(10.0, 60.0, length = n + 1) # linearly space wages
 
 using StatsPlots
-plt = plot(w, pdf.(dist, support(dist)), xlabel = "wages",
+plt = plot(w, p, xlabel = "wages",
            ylabel = "probabilities", legend = false)
 ```
 
 We can explore taking expectations over this distribution
 
 ```{code-cell} julia
-E = expectation(dist) # expectation operator
-
-# exploring the properties of the operator
+# exploring expectations via direct dot products
 wage(i) = w[i + 1] # +1 to map from support of 0
-E_w = E(wage)
-E_w_2 = E(i -> wage(i)^2) - E_w^2 # variance
+w = wage.(support(dist))
+E_w = dot(p, w)
+E_w_2 = dot(p, w .^ 2) - E_w^2 # variance
 @show E_w, E_w_2
 
-# use operator with left-multiply
-@show E * w # the `w` are values assigned for the discrete states
-@show dot(pdf.(dist, support(dist)), w); # identical calculation
+# could also use sum with elementwise product
+@show sum(p .* w);
 ```
 
 To implement our algorithm, let's have a look at the sequence of approximate value functions that
@@ -351,8 +350,8 @@ num_plots = 6
 
 # Operator
 # Broadcasts over the w, fixes the v
-T(v) = max.(w / (1 - beta), c + beta * E * v)
-# alternatively, T(v) = [max(wval/(1 - beta), c + beta * E*v) for wval in w]
+T(v) = max.(w / (1 - beta), c + beta * dot(v, p))
+# alternatively, T(v) = [max(w_val/(1 - beta), c + beta * dot(v, p)) for w_val in w]
 
 # fill in  matrix of vs
 vs = zeros(n + 1, 6) # data to fill
@@ -376,7 +375,7 @@ function compute_reservation_wage_direct(params;
     (; c, beta, w) = params
 
     # create a closure for the T operator
-    T(v) = max.(w / (1 - beta), c + beta * E * v) # (5) fixing the parameter values
+    T(v) = max.(w / (1 - beta), c + beta * dot(v, p)) # (5) fixing the parameter values
 
     v = copy(v_iv) # copy to prevent v_iv modification
     v_next = similar(v)
@@ -389,7 +388,7 @@ function compute_reservation_wage_direct(params;
         v .= v_next  # copy contents into v
     end
     # now compute the reservation wage
-    return (1 - beta) * (c + beta * E * v) # (2)
+    return (1 - beta) * (c + beta * dot(v, p)) # (2)
 end
 ```
 
@@ -413,12 +412,12 @@ In this case, we can use the `fixedpoint` algorithm discussed in {doc}`our Julia
 function compute_reservation_wage(params; v_iv = collect(w ./ (1 - beta)),
                                   iterations = 500, ftol = 1e-6, m = 1)
     (; c, beta, w) = params
-    T(v) = max.(w / (1 - beta), c + beta * E * v) # (5) fixing the parameter values
+    T(v) = max.(w / (1 - beta), c + beta * dot(v, p)) # (5) fixing the parameter values
 
     sol = fixedpoint(T, v_iv; iterations, ftol, m) # (5)
     sol.f_converged || error("Failed to converge")
     v_star = sol.zero
-    return (1 - beta) * (c + beta * E * v_star) # (3)
+    return (1 - beta) * (c + beta * dot(v_star, p)) # (3)
 end
 ```
 
@@ -578,9 +577,9 @@ The big difference here, however, is that we're iterating on a single number, ra
 Here's an implementation:
 
 ```{code-cell} julia
-function compute_reservation_wage_psi(c, beta; psi_iv = E * w ./ (1 - beta),
+function compute_reservation_wage_psi(c, beta; psi_iv = dot(w, p) / (1 - beta),
                                       iterations = 500, ftol = 1e-5, m = 1)
-    T_psi(psi) = [c + beta * E * max.((w ./ (1 - beta)), psi[1])] # (7)
+    T_psi(psi) = [c + beta * dot(max.((w ./ (1 - beta)), psi[1]), p)] # (7)
     # using vectors since fixedpoint doesn't support scalar
     sol = fixedpoint(T_psi, [psi_iv]; iterations, ftol, m)
     sol.f_converged || error("Failed to converge")
@@ -595,9 +594,9 @@ You can use this code to solve the exercise below.
 Another option is to solve for the root of the  $T_{\psi}(\psi) - \psi$ equation
 
 ```{code-cell} julia
-function compute_reservation_wage_psi2(c, beta; psi_iv = E * w ./ (1 - beta),
+function compute_reservation_wage_psi2(c, beta; psi_iv = dot(w, p) / (1 - beta),
                                        maxiters = 500, rtol = 1e-5)
-    root_psi(psi) = c + beta * E * max.((w ./ (1 - beta)), psi) - psi # (7)
+    root_psi(psi) = c + beta * dot(max.((w ./ (1 - beta)), psi), p) - psi # (7)
     psi_star = find_zero(root_psi, psi_iv; maxiters, rtol)
     return (1 - beta) * psi_star # (2)
 end

lectures/dynamic_programming/mccall_model_with_separation.md

@@ -47,7 +47,7 @@ Once separation enters the picture, the agent comes to view
 tags: [hide-output]
 ---
 using LinearAlgebra, Statistics
-using Distributions, Expectations, LaTeXStrings,NLsolve, Plots
+using Distributions, LaTeXStrings, NLsolve, Plots
 ```
 
 ## The Model
@@ -233,12 +233,12 @@ using Test
 ```
 
 ```{code-cell} julia
-using Distributions, LinearAlgebra, Expectations, NLsolve, Plots
+using Distributions, LinearAlgebra, NLsolve, Plots
 
 function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)),
                             tol = 1e-5,
                             iter = 2_000)
-    (; alpha, beta, sigma, c, gamma, w, dist, u) = mcm
+    (; alpha, beta, sigma, c, gamma, w, dist, u, p) = mcm
 
     # parameter validation
     @assert c > 0.0
@@ -247,14 +247,13 @@ function solve_mccall_model(mcm; U_iv = 1.0, V_iv = ones(length(mcm.w)),
     # necessary objects
     u_w = mcm.u.(w, sigma)
     u_c = mcm.u(c, sigma)
-    E = expectation(dist) # expectation operator for wage distribution
 
     # Bellman operator T. Fixed point is x* s.t. T(x*) = x*
     function T(x)
         V = x[1:(end - 1)]
         U = x[end]
         [u_w + beta * ((1 - alpha) * V .+ alpha * U);
-         u_c + beta * (1 - gamma) * U + beta * gamma * E * max.(U, V)]
+         u_c + beta * (1 - gamma) * U + beta * gamma * dot(max.(U, V), p)]
     end
 
     # value function iteration
@@ -293,7 +292,8 @@ function McCallModel(; alpha = 0.2,
                      u = (c, sigma) -> (c^(1 - sigma) - 1) / (1 - sigma),
                      w = range(10, 20, length = 60), # wage values
                      dist = BetaBinomial(59, 600, 400)) # distribution over wage values
-    return (; alpha, beta, gamma, c, sigma, u, w, dist)
+    p = pdf.(dist, support(dist))
+    return (; alpha, beta, gamma, c, sigma, u, w, dist, p)
 end
 ```
 

lectures/introduction_dynamics/kalman.md

@@ -643,44 +643,6 @@ tags: [remove-cell]
 end
 ```
 
-### Exercise 2
-
-```{code-cell} julia
-using Random, Expectations
-Random.seed!(43)  # reproducible results
-epsilon = 0.1
-kalman = Kalman(A, G, Q, R)
-set_state!(kalman, x_hat_0, Sigma_0)
-nodes, weights = qnwlege(21, theta - epsilon, theta + epsilon)
-
-T = 600
-z = zeros(T)
-for t in 1:T
-    # record the current predicted mean and variance, and plot their densities
-    m, v = kalman.cur_x_hat, kalman.cur_sigma
-    dist = Truncated(Normal(m, sqrt(v)), theta - 30 * epsilon, theta + 30 * epsilon) # define on compact interval, so we can use gridded expectation
-    E = expectation(dist, nodes) # nodes ∈ [theta-epsilon, theta+epsilon]
-    integral = E(x -> 1) # just take the pdf integral
-    z[t] = 1.0 - integral
-    # generate the noisy signal and update the Kalman filter
-    update!(kalman, theta + randn())
-end
-
-plot(1:T, z, fillrange = 0, color = :blue, fillalpha = 0.2, grid = false,
-     xlims = (0, T),
-     legend = false)
-```
-
-```{code-cell} julia
----
-tags: [remove-cell]
----
-@testset "Solution 2 Tests" begin
-    # #test z[4] ≈ 0.9310333042533682
-    @test T == 600
-end
-```
-
 ### Exercise 3
 
 ```{code-cell} julia

lectures/more_julia/general_packages.md

@@ -44,7 +44,7 @@ Also see {doc}`data and statistical packages <../more_julia/data_statistical_pac
 tags: [hide-output]
 ---
 using LinearAlgebra, Statistics
-using QuadGK, FastGaussQuadrature, Distributions, Expectations
+using QuadGK, FastGaussQuadrature, Distributions
 using Interpolations, Plots,  ProgressMeter
 ```
 
@@ -82,26 +82,6 @@ f(x) = x^2
 
 The only problem with the `FastGaussQuadrature` package is that you will need to deal with affine transformations to the non-default domains yourself.
 
-### Expectations
-
-If the calculations of the numerical integral is simply for calculating mathematical expectations of a particular distribution, then [Expectations.jl](https://github.com/QuantEcon/Expectations.jl) provides a convenient interface.
-
-Under the hood, it is finding the appropriate Gaussian quadrature scheme for the distribution using `FastGaussQuadrature`.
-
-```{code-cell} julia
-using Distributions, Expectations
-dist = Normal()
-E = expectation(dist)
-f(x) = x
-@show E(f) #i.e. identity
-
-# Or using as a linear operator
-f(x) = x^2
-x = nodes(E)
-w = weights(E)
-E * f.(x) == f.(x) ⋅ w
-```
-
 ## Interpolation
 
 In economics we often wish to interpolate discrete data (i.e., build continuous functions that join discrete sequences of points).