Merge 0d03157 into 294a551

.travis.yml

@@ -1,7 +1,6 @@
 language: julia
 sudo: false
 julia:
-    - 0.4
     - 0.5
     - nightly
 matrix:

REQUIRE

@@ -1,4 +1,4 @@
-julia 0.4
+julia 0.5
 DSP
 Distributions 0.8.9
 LightGraphs

src/QuantEcon.jl

@@ -11,8 +11,12 @@ using DSP: TFFilter, freqz
 using Primes: primes
 using Compat: view, @compat
 
+@static if isdefined(Base, :Iterators)
+    using Base.Iterators: cycle, take
+end
+
 # useful types
-typealias ScalarOrArray{T} Union{T, Array{T}}
+@compat ScalarOrArray{T} = Union{T,Array{T}}
 
 
 export

src/arma.jl

@@ -119,7 +119,7 @@ function spectral_density(arma::ARMA; res=1200, two_pi::Bool=true)
     w = linspace(0, wmax, res)
     tf = TFFilter(reverse(arma.ma_poly), reverse(arma.ar_poly))
     h = freqz(tf, w)
-    spect = arma.sigma^2 * abs(h).^2
+    spect = arma.sigma.^2 .* abs.(h).^2
     return w, spect
 end
 
@@ -166,7 +166,7 @@ function impulse_response(arma::ARMA; impulse_length=30)
     # == Pad theta with zeros at the end == #
     theta = [arma.theta; zeros(impulse_length - arma.q)]
     psi_zero = 1.0
-    psi = Array(Float64, impulse_length)
+    psi = Array{Float64}(impulse_length)
     for j = 1:impulse_length
         psi[j] = theta[j]
         for i = 1:min(j, arma.p)
@@ -197,7 +197,7 @@ function simulation(arma::ARMA; ts_length=90, impulse_length=30)
     T = ts_length
     psi = impulse_response(arma, impulse_length=impulse_length)
     epsilon = arma.sigma * randn(T + J)
-    X = Array(Float64, T)
+    X = Array{Float64}(T)
     for t=1:T
         X[t] = dot(epsilon[t:J+t-1], psi)
     end

src/compute_fp.jl

@@ -64,7 +64,7 @@ function compute_fixed_point{TV}(T::Function,
     while iterate < max_iter && err > err_tol
         new_v = T(v)::TV
         iterate += 1
-        err = Base.maxabs(new_v - v)
+        err = Base.maximum(abs, new_v - v)
         if verbose == 2
             if iterate % print_skip == 0
                 println("Compute iterate $iterate with error $err")

src/ecdf.jl

@@ -21,15 +21,14 @@ type ECDF
     observations::Vector
 end
 
-ecdf(e::ECDF, x::Real) = mean(e.observations .<= x)
-ecdf(e::ECDF, x::Array) = map(i->ecdf(e, i), x)
-
 """
 Evaluate the empirical cdf at one or more points
 
 ##### Arguments
 
-- `e::ECDF`: The `ECDF` instance
 - `x::Union{Real, Array}`: The point(s) at which to evaluate the ECDF
 """
-ecdf
+(e::ECDF)(x::Real) = mean(e.observations .<= x)
+(e::ECDF)(x::AbstractArray) = e.(x)
+
+@deprecate ecdf(e::ECDF, x) e(x)

src/estspec.jl

@@ -73,7 +73,7 @@ end
 
 function periodogram(x::Vector)
     n = length(x)
-    I_w = abs(fft(x)).^2 ./ n
+    I_w = abs.(fft(x)).^2 ./ n
     w = 2pi * (0:n-1) ./ n  # Fourier frequencies
 
     # int rounds to nearest integer. We want to round up or take 1/2 + 1 to
@@ -148,7 +148,7 @@ function ar_periodogram(x::Array, window::AbstractString="hanning", window_len::
     w, I_w = periodogram(e_hat, window, window_len)
 
     # recolor and return
-    I_w = I_w ./ abs(1 - phi .* exp(im.*w)).^2
+    I_w = I_w ./ abs.(1 .- phi .* exp.(im.*w)).^2
 
     return w, I_w
 end

src/lqcontrol.jl

@@ -95,14 +95,14 @@ symmetric and nonnegative definite
 Must be symmetric and nonnegative definite
 - `A::ScalarOrArray` : n x n coefficient on state in state transition
 - `B::ScalarOrArray` : n x k coefficient on control in state transition
-- `;C::ScalarOrArray(zeros(size(R, 1)))` : n x j coefficient on random shock in
+- `;C::ScalarOrArray{zeros(size(R}(1)))` : n x j coefficient on random shock in
 state transition
-- `;N::ScalarOrArray(zeros(size(B,1), size(A, 2)))` : k x n cross product in
+- `;N::ScalarOrArray{zeros(size(B,1)}(size(A, 2)))` : k x n cross product in
 payoff equation
 - `;bet::Real(1.0)` : Discount factor in [0, 1]
 - `capT::Union{Int, Void}(Void)` : Terminal period in finite horizon
 problem
-- `rf::ScalarOrArray(fill(NaN, size(R)...))` : n x n terminal payoff in finite
+- `rf::ScalarOrArray{fill(NaN}(size(R)...))` : n x n terminal payoff in finite
 horizon problem. Must be symmetric and nonnegative definite.
 
 """
@@ -232,8 +232,8 @@ Private method implementing `compute_sequence` when state is a scalar
 function _compute_sequence{T}(lq::LQ, x0::T, policies)
     capT = length(policies)
 
-    x_path = Array(T, capT+1)
-    u_path = Array(T, capT)
+    x_path = Array{T}(capT+1)
+    u_path = Array{T}(capT)
 
     x_path[1] = x0
     u_path[1] = -(first(policies)*x0)
@@ -259,8 +259,8 @@ function _compute_sequence{T}(lq::LQ, x0::Vector{T}, policies)
 
     A, B, C = lq.A, reshape(lq.B, n, k), reshape(lq.C, n, j)
 
-    x_path = Array(T, n, capT+1)
-    u_path = Array(T, k, capT)
+    x_path = Array{T}(n, capT+1)
+    u_path = Array{T}(k, capT)
     w_path = C*randn(j, capT+1)
 
     x_path[:, 1] = x0
@@ -305,7 +305,7 @@ function compute_sequence(lq::LQ, x0::ScalarOrArray, ts_length::Integer=100)
         policies = fill(lq.F, ts_length)
     else
         capT = min(ts_length, lq.capT)
-        policies = Array(typeof(lq.F), capT)
+        policies = Array{typeof(lq.F)}(capT)
         for t = capT:-1:1
             update_values!(lq)
             policies[t] = lq.F

src/lqnash.jl

@@ -100,7 +100,7 @@ function nnash(a, b1, b2, r1, r2, q1, q2, s1, s2, w1, w2, m1, m2;
         p1 = (a2'*p1*a2) .+ r1 .+ (f2'*s1*f2) .- (a2'*p1*b1 .+ w1 .- f2'*m1)*f1
         p2 = (a1'*p2*a1) .+ r2 .+ (f1'*s2*f1) .- (a1'*p2*b2 .+ w2 .- f1'*m2)*f2
 
-        dd = maximum(abs(f10 .- f1) + abs(f20 .- f2))
+        dd = maximum(abs.(f10 .- f1) + abs.(f20 .- f2))
         its = its + 1
         if its > max_iter
             error("Reached max iterations, no convergence")

src/lss.jl

@@ -111,7 +111,7 @@ end
 
 
 function simulate(lss::LSS, ts_length=100)
-    x = Array(Float64, lss.n, ts_length)
+    x = Array{Float64}(lss.n, ts_length)
     x[:, 1] = rand(lss.dist)
     w = randn(lss.m, ts_length - 1)
     for t=1:ts_length-1
@@ -140,7 +140,7 @@ Simulate num_reps observations of x_T and y_T given x_0 ~ N(mu_0, Sigma_0).
 
 """
 function replicate(lss::LSS, t::Integer, num_reps::Integer=100)
-    x = Array(Float64, lss.n, num_reps)
+    x = Array{Float64}(lss.n, num_reps)
     for j=1:num_reps
         x_t, _ = simulate(lss, t+1)
         x[:, j] = x_t[:, end]
@@ -153,8 +153,27 @@ end
 replicate(lss::LSS; t::Integer=10, num_reps::Integer=100) =
     replicate(lss, t, num_reps)
 
+immutable LSSMoments
+    lss::LSS
+end
+
+Base.start(L::LSSMoments) = (copy(L.lss.mu_0), copy(L.lss.Sigma_0))
+Base.done(L::LSSMoments, _) = false
+function Base.next(L::LSSMoments, moms)
+    A, C, G = L.lss.A, L.lss.C, L.lss.G
+    mu_x, Sigma_x = moms
+
+    mu_y, Sigma_y = G * mu_x, G * Sigma_x * G'
+
+    # Update moments of x
+    mu_x2 = A * mu_x
+    Sigma_x2 = A * Sigma_x * A' + C * C'
+
+    (mu_x, mu_y, Sigma_x, Sigma_y), (mu_x2, Sigma_x2)
+end
+
 """
-Create a generator to calculate the population mean and
+Create an iterator to calculate the population mean and
 variance-convariance matrix for both x_t and y_t, starting at
 the initial condition (self.mu_0, self.Sigma_0).  Each iteration
 produces a 4-tuple of items (mu_x, mu_y, Sigma_x, Sigma_y) for
@@ -165,18 +184,8 @@ the next period.
 - `lss::LSS` An instance of the Gaussian linear state space model
 
 """
-function moment_sequence(lss::LSS)
-    A, C, G = lss.A, lss.C, lss.G
-    mu_x, Sigma_x = copy(lss.mu_0), copy(lss.Sigma_0)
-    while true
-        mu_y, Sigma_y = G * mu_x, G * Sigma_x * G'
-        produce((mu_x, mu_y, Sigma_x, Sigma_y))
-
-        # Update moments of x
-        mu_x = A * mu_x
-        Sigma_x = A * Sigma_x * A' + C * C'
-    end
-end
+moment_sequence(lss::LSS) = LSSMoments(lss)
+
 
 """
 Compute the moments of the stationary distributions of x_t and
@@ -199,22 +208,22 @@ initial conditions lss.mu_0 and lss.Sigma_0
 """
 function stationary_distributions(lss::LSS; max_iter=200, tol=1e-5)
     # Initialize iteration
-    m = @task moment_sequence(lss)
-    mu_x, mu_y, Sigma_x, Sigma_y = consume(m)
+    m = moment_sequence(lss)
+    mu_x, mu_y, Sigma_x, Sigma_y = first(m)
 
     i = 0
-    err = tol + 1.
-
-    while err > tol
-        if i > max_iter
-            error("Convergence failed after $i iterations")
-        else
-            i += 1
-            mu_x1, mu_y, Sigma_x1, Sigma_y = consume(m)
-            err_mu = Base.maxabs(mu_x1 - mu_x)
-            err_Sigma = Base.maxabs(Sigma_x1 - Sigma_x)
-            err = max(err_Sigma, err_mu)
-            mu_x, Sigma_x = mu_x1, Sigma_x1
+    err = tol + 1.0
+
+    for (mu_x1, mu_y, Sigma_x1, Sigma_y) in m
+        i > max_iter && error("Convergence failed after $i iterations")
+        i += 1
+        err_mu = maximum(abs, mu_x1 - mu_x)
+        err_Sigma = maximum(abs, Sigma_x1 - Sigma_x)
+        err = max(err_Sigma, err_mu)
+        mu_x, Sigma_x = mu_x1, Sigma_x1
+
+        if err < tol && i > 1
+            break
         end
     end
 

src/markov/ddp.jl

@@ -115,7 +115,7 @@ type DiscreteDP{T<:Real,NQ,NR,Tbeta<:Real,Tind}
         end
 
         if _has_sorted_sa_indices(s_indices, a_indices)
-            a_indptr = Array(Int64, num_states+1)
+            a_indptr = Array{Int64}(num_states+1)
             _a_indices = copy(a_indices)
             _generate_a_indptr!(num_states, s_indices, a_indptr)
         else
@@ -539,7 +539,7 @@ end
 
 # TODO: express it in a similar way as above to exploit Julia's column major order
 function RQ_sigma{T<:Integer}(ddp::DDPsa, sigma::Array{T})
-    sigma_indices = Array(T, num_states(ddp))
+    sigma_indices = Array{T}(num_states(ddp))
     _find_indices!(get(ddp.a_indices), get(ddp.a_indptr), sigma, sigma_indices)
     R_sigma = ddp.R[sigma_indices]
     Q_sigma = ddp.Q[sigma_indices, :]
@@ -602,7 +602,7 @@ end
 
 function s_wise_max(ddp::DDPsa, vals::Vector)
     s_wise_max!(get(ddp.a_indices), get(ddp.a_indptr), vals,
-                 Array(Float64, num_states(ddp)))
+                 Array{Float64}(num_states(ddp)))
 end
 
 s_wise_max!(ddp::DDPsa, vals::Vector, out::Vector, out_argmax::Vector) =
@@ -751,7 +751,7 @@ function _solve!(ddp::DiscreteDP, ddpr::DPSolveResult{VFI}, max_iter::Integer,
         bellman_operator!(ddp, ddpr)
 
         # compute error and update the v inside ddpr
-        err = maxabs(ddpr.Tv .- ddpr.v)
+        err = maximum(abs, ddpr.Tv .- ddpr.v)
         copy!(ddpr.v, ddpr.Tv)
         ddpr.num_iter += 1
 

src/markov/mc_tools.jl

@@ -13,7 +13,7 @@ http://quant-econ.net/jl/finite_markov.html
 import LightGraphs: DiGraph, period, attracting_components,
                     strongly_connected_components, is_strongly_connected
 
-@inline check_stochastic_matrix(P) = maxabs(sum(P, 2) - 1) < 1e-15 ? true : false
+@inline check_stochastic_matrix(P) = maximum(abs, sum(P, 2) - 1) < 1e-15 ? true : false
 
 """
 Finite-state discrete-time Markov chain.
@@ -240,7 +240,7 @@ for (S, ex_T, ex_gth) in ((Real, :(T), :(gth_solve!)),
         n = n_states(mc)
         rec_classes = recurrent_classes(mc)
         T1 = $ex_T
-        stationary_dists = Array(Vector{T1}, length(rec_classes))
+        stationary_dists = Array{Vector{T1}}(length(rec_classes))
 
         for (i, rec_class) in enumerate(rec_classes)
             dist = zeros(T1, n)
@@ -361,7 +361,7 @@ ts_length
 """
 function simulate(mc::MarkovChain, ts_length::Int;
                   init::Int=rand(1:n_states(mc)))
-    X = Array(eltype(mc), ts_length)
+    X = Array{eltype(mc)}(ts_length)
     simulate!(X, mc; init=init)
 end
 
@@ -416,7 +416,7 @@ ts_length
 """
 function simulate_indices(mc::MarkovChain, ts_length::Int;
                           init::Int=rand(1:n_states(mc)))
-    X = Array(Int, ts_length)
+    X = Array{Int}(ts_length)
     simulate_indices!(X, mc; init=init)
 end
 

src/markov/random_mc.jl

@@ -215,7 +215,7 @@ function random_probvec(k::Integer, m::Integer)
     k == 1 && return ones((k, m))
 
     # if k >= 2
-    x = Array(Float64, (k, m))
+    x = Array{Float64}((k, m))
 
     r = rand(k-1, m)
     x[1:end-1, :] = sort(r, 1)

src/matrix_eqn.jl

@@ -46,7 +46,7 @@ function solve_discrete_lyapunov(A::ScalarOrArray,
         alpha1 = alpha0*alpha0
         gamma1 = gamma0 + alpha0*gamma0*alpha0'
 
-        diff = maximum(abs(gamma1 - gamma0))
+        diff = maximum(abs, gamma1 - gamma0)
         alpha0 = alpha1
         gamma0 = gamma1
 
@@ -159,7 +159,7 @@ function solve_discrete_riccati(A::ScalarOrArray, B::ScalarOrArray,
         G1 = G0 + A0 * G0 * ((I + H0 * G0)\A0')
         H1 = H0 + A0' * ((I + H0*G0)\(H0*A0))
 
-        dist = Base.maxabs(H1 - H0)
+        dist = Base.maximum(abs, H1 - H0)
         A0 = A1
         G0 = G1
         H0 = H1

src/optimization.jl

@@ -19,7 +19,7 @@ function golden_method(f::Function, a::AbstractVector, b::AbstractVector;
         f1[i] = f2[i]
         d *= α2
         x2 = x1 + s.*(i- ~i).*d
-        s = sign(x2 - x1)
+        s = sign.(x2 .- x1)
         f2 = f(x2)
     end