Ps2

.github/workflows/docs.yml

@@ -0,0 +1,22 @@
+name: Documenter
+on:
+  push:
+    branches: [main]
+    tags: [v*]
+  pull_request:
+
+jobs:
+  Documenter:
+    permissions:
+      contents: write
+    name: Documentation
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+      - uses: julia-actions/cache@v1
+        with:
+          cache-compiled: "true"
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-docdeploy@v1
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

.github/workflows/tests.yml

@@ -0,0 +1,30 @@
+name: Tests
+on:
+  push:
+    branches: ['main']
+    tags: '*'
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: ['1']
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/cache@v1
+        with:
+          cache-compiled: "true"
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
+

.gitignore

@@ -1,3 +1,6 @@
+# emacs backups
+*~
+
 # Files generated by invoking Julia with --code-coverage
 *.jl.cov
 *.jl.*.cov

Project.toml

@@ -1,7 +1,15 @@
+name = "ECON622_BLP"
+uuid = "04c43c4b-146f-4579-99e2-63279d378c20"
+authors = ["Paul Schrimpf <[email protected]>"]
+version = "0.0.1"
+
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Sobol = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
 SparseGrids = "bafbe729-afc6-5148-bb4f-226bf3d46895"

docs/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

docs/make.jl

@@ -0,0 +1,18 @@
+using Documenter, ECON622_BLP
+makedocs(;
+         modules=[ECON622_BLP],
+         format=Documenter.HTML(),
+         pages=[
+           "Home" => "index.md",
+           "Function Reference" => "functions.md"
+         ],
+         repo=Remotes.GitHub("UBCECON","ECON622_BLP.jl"),
+         sitename="ECON622_BLP.jl",
+         #doctest=false,
+         authors="Paul Schrimpf <[email protected]>"
+         #assets=String[],
+         #plugins=[bib]
+)
+
+
+deploydocs(repo="github.com/UBCECON/ECON622_BLP.jl.git", branch="gh-pages")

docs/src/functions.md

@@ -0,0 +1,10 @@
+# Function Reference
+
+```@autodocs
+Modules = [ECON622_BLP]
+```
+
+## Index
+
+```@index
+```

docs/src/index.md

@@ -0,0 +1,10 @@
+# ECON622_BLP.jl
+
+Random coefficients demand models in the style of Berry, Levinsohn, and Pakes (1995). Created as an example for teaching ECON622.
+
+```@contents
+```
+
+# Model
+
+# Another Section

src/ECON622_BLP.jl

@@ -0,0 +1,123 @@
+module ECON622_BLP
+
+import NLsolve
+
+export share, delta, Integrate, share!
+
+include("integrate.jl")
+
+@doc raw"""
+    share(δ, Σ, dFν, x)
+
+Computes shares in random coefficient logit with mean tastes `δ`, observed characteristics `x`, unobserved taste distribution `dFν`, and taste covariances `Σ`. 
+Assumes there is an outside option with u=0. The outside option has share `1-sum(s)`
+
+# Arguments
+
+- `δ` vector of length `J`
+- `Σ` `K` by `K` matrix
+- `x` `J` by `K` array
+- `∫` AbstractIntegrator for integrating over distribution of `ν`
+
+# Returns
+
+- vector of length `J` consisting of $s_1$, ..., $s_J$
+"""
+function share(δ, Σ, x, ∫::Integrate.AbstractIntegrator)
+  J,K = size(x)
+  (length(δ) == J) || error("length(δ)=$(length(δ)) != size(x,1)=$J")
+  (K,K) == size(Σ) || error("size(x,2)=$K != size(Σ)=$(size(Σ))")
+  function shareν(ν)
+    s = δ + x*Σ*ν
+    smax=max(0,maximum(s))    
+    return(exp.(s .- smax)./(sum(exp,s .- smax) + exp(-smax)))
+  end
+  return(∫(shareν))
+end
+
+function delta(s, Σ, x, ∫)
+  δ = similar(s)
+  delta!(δ, s, Σ, x, ∫)
+  return(δ)
+end
+
+# to define custom reverse rule for enzyme
+function delta!(δ, s, Σ, x, ∫)
+  eqns! = let xi=x, Σi=Σ, ∫i=∫ 
+    function(F,δ) 
+        F .= s - share(δ,Σi,xi,∫i)
+    end
+  end
+  δ0 = log.(s) .- log(1-sum(s))
+  sol=NLsolve.nlsolve(eqns!, δ0, autodiff=:forward, method=:trust_region)
+  if (sol.residual_norm > 1e-4)
+    @warn "Possible problem in delta!(...)\n".*"$sol"
+  end
+  δ .= sol.zero
+  return nothing
+end 
+
+
+import EnzymeCore 
+import EnzymeCore.EnzymeRules
+import EnzymeCore: Const, Duplicated, BatchDuplicated, Active, BatchDuplicatedNoNeed
+import EnzymeCore.EnzymeRules: ConfigWidth, AugmentedReturn, needs_primal, needs_shadow
+import Enzyme # needed for Enzyme.jacobian
+import LinearAlgebra: lu
+
+function EnzymeRules.forward(f::Const{typeof(delta)}, ::Type{<:Duplicated}, 
+  s::Duplicated, Σ::Const, x::Const, ∫::Const)
+  δ  = f.val(s.val, Σ.val, x.val, ∫.val)
+  sf(d) = share(d, Σ.val, x.val, ∫.val)
+  J = Enzyme.jacobian(Enzyme.Forward,sf, δ)
+  dδ = J \ s.dval    
+  return Duplicated(δ, dδ)
+end
+
+function EnzymeRules.forward(f::Const{typeof(delta)}, RT::Type{<:BatchDuplicated}, 
+                             s::BatchDuplicated, Σ::Const, x::Const, ∫::Const)
+  δ  = f.val(s.val, Σ.val, x.val, ∫.val)
+  sf(d) = share(d, Σ.val, x.val, ∫.val)
+  J = lu(Enzyme.jacobian(Enzyme.Forward,sf, δ))
+  out = BatchDuplicated(δ, map(ds->J \ ds, s.dval))
+  return out
+end
+
+function EnzymeRules.forward(f::Const{typeof(delta)}, RT::Type{<:BatchDuplicatedNoNeed}, 
+                             s::BatchDuplicated{T,N}, Σ::Const, x::Const, ∫::Const) where {T, N}
+  δ  = f.val(s.val, Σ.val, x.val, ∫.val)
+  sf(d) = share(d, Σ.val, x.val, ∫.val)
+  J = lu(Enzyme.jacobian(Enzyme.Forward,sf, δ))
+  out = NTuple{N,T}(J \ ds for ds in s.dval)
+  return out
+end
+
+
+function EnzymeRules.augmented_primal(config::ConfigWidth{W}, f::Const{typeof(delta!)}, RT::Type{<:Const}, 
+                                      δ::Duplicated, s::Duplicated, Σ::Const, x::Const, ∫) where W
+  f.val(δ.val,s.val, Σ.val, x.val, ∫.val)
+  tape = (δ = copy(δ.val),s=copy(s.val),Σ=copy(Σ.val))
+  @assert !needs_primal(config) && !needs_shadow(config)
+  return AugmentedReturn(nothing , nothing, tape)
+end
+
+function EnzymeRules.reverse(config::ConfigWidth{W}, f::Const{typeof(delta!)}, dret::Type{<:Const}, tape, 
+                             δ::Duplicated, s::Duplicated, Σ::Const, x::Const, ∫) where W
+  sf(d) =share(d, tape.Σ, x.val, ∫.val)
+  J = lu(Enzyme.jacobian(Enzyme.Forward,sf, tape.δ))
+  s.dval .+= J \ δ.dval 
+  return(nothing, nothing, nothing,nothing,nothing)
+end
+
+function broken!(x)
+  out = 0
+  for i ∈ eachindex(x)
+    out+=exp(x[i]) 
+  end
+  x[1] = x[length(x)+1]
+  return(out)
+end
+
+
+
+end

src/integrate.jl

@@ -2,36 +2,57 @@ module Integrate
 
 using Distributions
 import Sobol: skip, SobolSeq
-import Base.Iterators: take, Repeated
+import Base.Iterators: take, Repeated, product, repeated
 import HCubature: hcubature
 import LinearAlgebra: cholesky
+using FastGaussQuadrature, LinearAlgebra
+using SparseGrids
+
 
 abstract type AbstractIntegrator end
 
-(∫::AbstractIntegrator)(f::Function) = sum(w*f(x) for (w,x) in zip(∫.w, ∫.x))
+(∫::AbstractIntegrator)(f::Function) = #sum(∫.w[i]*f(∫.x[i]) for i in eachindex(∫.x))
+    sum(w*f(x) for (x,w) in zip(∫.x, ∫.w))
 
 struct FixedNodeIntegrator{Tx,Tw} <: AbstractIntegrator
     x::Tx
     w::Tw
 end
 
-MonteCarloIntegrator(distribution::Distribution, ndraw=100)=FixedNodeIntegrator([rand(distribution) for i=1:ndraw], Repeated(1/ndraw))
+MonteCarloIntegrator(distribution::Distribution, ndraw=100)=
+  FixedNodeIntegrator([rand(distribution) for i=1:ndraw], fill(1/ndraw, ndraw))
 
 function QuasiMonteCarloIntegrator(distribution::UnivariateDistribution, ndraws=100)
     ss = skip(SobolSeq(1), ndraw)
     x = [quantile(distribution, x) for x in take(ss,ndraw)]
-    w = Repeated(1/ndraw)
+    w = fill(1/ndraw, ndraw) #repeated(1/ndraw, ndraw)
     FixedNodeIntegrator(x,w)
 end 
 
 function QuasiMonteCarloIntegrator(distribution::AbstractMvNormal, ndraw=100)
     ss = skip(SobolSeq(length(distribution)), ndraw)
     L = cholesky(distribution.Σ).L
     x = [L*quantile.(Normal(), x) for x in take(ss,ndraw)]
-    w = Repeated(1/ndraw)
+    w = fill(1/ndraw, ndraw) #repeated(1/ndraw, ndraw)
     FixedNodeIntegrator(x,w)
 end 
 
+function  QuadratureIntegrator(dx::AbstractMvNormal; ndraw=100)
+  n = Int(ceil(ndraw^(1/length(dx))))
+  x, w = gausshermite(n)
+  L = cholesky(dx.Σ).L
+  x = vec([√2*L*vcat(xs...) + dx.μ for xs in product(repeated(x,length(dx))...)])
+  w = vec([prod(ws)/π^(length(dx)/2) for ws in product(repeated(w, length(dx))...)])
+  FixedNodeIntegrator(x,w)
+end
+
+function SparseGridIntegrator(dx::AbstractMvNormal; order=5)
+  X, W = sparsegrid(length(dx), order, gausshermite, sym=true)
+  L = cholesky(dx.Σ).L
+  X = [√2*L*x + dx.μ for x ∈ X]
+  W /= π^(length(dx)/2)
+  FixedNodeIntegrator(X,W)
+end
 
 struct AdaptiveIntegrator{FE,FT,FJ,A,L} <: AbstractIntegrator
     eval::FE
@@ -41,8 +62,6 @@ struct AdaptiveIntegrator{FE,FT,FJ,A,L} <: AbstractIntegrator
     limits::L
 end
 
-(∫::AdaptiveIntegrator)(f::Function) = ∫.eval(t->f(∫.transform(t))*∫.detJ(t), ∫.limits...; ∫.args...)[1]
-
 function AdaptiveIntegrator(dist::AbstractMvNormal; eval=hcubature, options=())
     D = length(dist)
     x(t) = t./(1 .- t.^2)
@@ -52,4 +71,6 @@ function AdaptiveIntegrator(dist::AbstractMvNormal; eval=hcubature, options=())
     AdaptiveIntegrator(hcubature,x,Dx,args, limits)
 end
 
+(∫::AdaptiveIntegrator)(f::Function) = ∫.eval(t->f(∫.transform(t))*∫.detJ(t), ∫.limits...; ∫.args...)[1]
+
 end

test/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+AbbreviatedStackTraces = "ac637c84-cc71-43bf-9c33-c1b4316be3d4"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"