Merge 67328c5 into 136ea7e

.travis.yml

@@ -1,13 +1,19 @@
 language: julia
+dist: trusty
+sudo: false
 os:
   - linux
 julia:
-  - 0.4
+  - 0.6
+  - nightly
 notifications:
   email: false
-sudo: false
-script:
-  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
-  - julia -e 'Pkg.clone(pwd()); Pkg.build("TimeModels"); Pkg.test("TimeModels", coverage=true)'
+addons:
+  apt:
+    packages:
+      - libnlopt0
+# script:
+  # - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
+  # - julia -e 'Pkg.clone(pwd()); Pkg.build("TimeModels"); Pkg.test("TimeModels", coverage=true)'
 after_success:
   - julia -e 'cd(Pkg.dir("TimeModels")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder())'

README.md

@@ -2,9 +2,9 @@ TimeModels.jl
 ============
 [![Build Status](https://travis-ci.org/JuliaStats/TimeModels.jl.svg?branch=master)](https://travis-ci.org/JuliaStats/TimeModels.jl)
 [![Coverage Status](https://coveralls.io/repos/JuliaStats/TimeModels.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaStats/TimeModels.jl?branch=master)
-[![TimeModels](http://pkg.julialang.org/badges/TimeModels_0.4.svg)](http://pkg.julialang.org/?pkg=TimeModels&ver=0.4)
+[![TimeModels](http://pkg.julialang.org/badges/TimeModels_0.6.svg)](http://pkg.julialang.org/?pkg=TimeModels&ver=0.6)
 
-##A Julia package for modeling time series. 
+##A Julia package for modeling time series.
 
 ![Kalman Demo](png/kalman.png)
 ![Experimental acf plot](png/acf_plot.png)
@@ -18,7 +18,7 @@ Generalized Autoregressive Conditional Heteroskedastic ([GARCH](http://en.wikipe
 * garchFit - estimates parameters of univariate normal GARCH process.
 * predict - make prediction using fitted object returned by garchFit
 * garchPkgTest - runs package test (compares model parameters with those obtained using R fGarch)
-* Jarque-Bera residuals test 
+* Jarque-Bera residuals test
 * Error analysis
 
 Analysis of model residuals - currently only Jarque-Bera Test implemented.
@@ -37,22 +37,22 @@ data - data vector
 #### returns:
 Structure containing details of the GARCH fit with the following fields:
 
-* data - orginal data  
-* params - vector of model parameters (omega,alpha,beta)  
-* llh - likelihood  
-* status - status of the solver  
-* converged - boolean convergence status, true if constraints are satisfied  
-* sigma - conditional volatility  
+* data - orginal data
+* params - vector of model parameters (omega,alpha,beta)
+* llh - likelihood
+* status - status of the solver
+* converged - boolean convergence status, true if constraints are satisfied
+* sigma - conditional volatility
 * hessian - Hessian matrix
 * secoef - standard errors
 * tval - t-statistics
-  
+
 ### predict
-make volatility prediction  
+make volatility prediction
 #### arguments:
-fit - fitted object returned by garchFit  
+fit - fitted object returned by garchFit
 #### returns:
-one-step-ahead volatility forecast  
+one-step-ahead volatility forecast
 
 ## Example
 

REQUIRE

@@ -1,4 +1,4 @@
-julia 0.4 0.5-
+julia 0.6
 Distributions
 StatsBase
 TimeSeries 0.5.11

src/GARCH.jl

@@ -39,11 +39,11 @@ function cdHessian(par,LLH)
   eps = 1e-4 * par
   n = length(par)
   H = zeros(n,n)
-  for(i = 1:n)
-    for(j = 1:n)
-      x1 = copy(par) 
+  for i = 1:n
+    for j = 1:n
+      x1 = copy(par)
       x1[i] += eps[i]
-      x1[j] += eps[j] 
+      x1[j] += eps[j]
       x2 = copy(par)
       x2[i] += eps[i]
       x2[j] -= eps[j]
@@ -61,7 +61,7 @@ end
 
 function garchLLH(rets::Vector,x::Vector)
   rets2   = rets.^2;
-  T = length(rets); 
+  T = length(rets);
   ht = zeros(T);
   omega,alpha,beta = x;
   ht[1] = sum(rets2)/T;
@@ -75,7 +75,7 @@ function predict(fit::GarchFit)
  omega, alpha, beta = fit.params;
  rets = fit.data
  rets2   = rets.^2;
- T = length(rets); 
+ T = length(rets);
  ht    = zeros(T);
  ht[1] = sum(rets2)/T;
  for i=2:T
@@ -87,7 +87,7 @@ end
 function garchFit(data::Vector)
   rets = data
   rets2   = rets.^2;
-  T = length(rets); 
+  T = length(rets);
   ht = zeros(T);
   function garchLike(x::Vector, grad::Vector)
     omega,alpha,beta = x;

src/TimeModels.jl

@@ -3,14 +3,15 @@ module TimeModels
 using Base.Dates
 using Distributions
 using StatsBase
-using TimeSeries 
+using TimeSeries
 using Optim
 using NLopt
+using Compat
 
 import Base: show
 
 export
-  # General state space model 
+  # General state space model
   StateSpaceModel, simulate,
 
   #Kalman
@@ -19,8 +20,8 @@ export
 
   # Parameter fitting
   ParametrizedMatrix,
-  parametrize_full, parametrize_diag, parametrize_none, 
-  ParametrizedSSM, SSMParameters, fit, 
+  parametrize_full, parametrize_diag, parametrize_none,
+  ParametrizedSSM, SSMParameters, fit,
 
   # ARIMA exports
   arima_statespace,
@@ -51,4 +52,4 @@ include("diagnostic_tests.jl")
 
 include("utilities.jl")
 
-end 
+end

src/kalman_filter.jl

@@ -26,7 +26,7 @@ function kalman_filter{T}(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=ze
     function kalman_recursions(y_i::Vector{T}, u_i::Vector{T},
                                   C_i::Matrix{T}, D_i::Matrix{T}, W_i::Matrix{T},
                                   x_pred_i::Vector{T}, P_pred_i::Matrix{T})
-        if !any(isnan(y_i))
+        if !any(isnan, y_i)
             innov =  y_i - C_i * x_pred_i - D_i * u_i
             S = C_i * P_pred_i * C_i' + W_i  # Innovation covariance
             K = P_pred_i * C_i' / S # Kalman gain
@@ -78,12 +78,12 @@ function loglikelihood{T}(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=ze
     y, u = y', u'
     n = size(y, 2)
 
-    @assert !any(isnan(u))
+    @assert !any(isnan, u)
     @assert size(y, 1) == model.ny
     @assert size(u, 1) == model.nu
     @assert size(u, 2) == n
 
-    y_notnan = !isnan(y)
+    y_notnan = (!).(isnan.(y))
     y = y .* y_notnan
 
     I0ny = zeros(model.ny, model.ny)

src/kalman_smooth.jl

@@ -45,14 +45,19 @@ function kalman_smooth{T}(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=ze
     u = u'
     n = size(y, 2)
 
-    @assert !any(isnan(u))
-    @assert size(y, 1) == model.ny
+    @assert !any(isnan, u)
+    # @assert size(y, 1) == model.ny
+    if size(y, 1) != model.ny
+        @show size(y)
+        @show model.ny
+        error("")
+    end
     @assert size(u, 1) == model.nu
     @assert size(u, 2) == n
 
     I0ny = zeros(model.ny, model.ny)
 
-    y_notnan = !isnan(y)
+    y_notnan = (!).(isnan.(y))
     y = y .* y_notnan
 
     x_pred_t, P_pred_t  = model.x1, model.P1
@@ -105,7 +110,7 @@ function kalman_smooth{T}(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=ze
         missing_obs = !all(y_notnan[:, t])
         if missing_obs
             ynnt = y_notnan[:, t]
-            I1, I2 = spdiagm(ynnt), spdiagm(!ynnt)
+            I1, I2 = spdiagm(ynnt), spdiagm(.!ynnt)
             Ct, Dut = I1 * Ct, I1 * Dut
             Wt = I1 * Wt * I1 + I2 * Wt * I2
         end #if

src/parameter_estimation.jl

@@ -17,7 +17,7 @@ function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
 
     y_orig = copy(y)
     y = y'
-    y_notnan = !isnan(y)
+    y_notnan = (!).(isnan.(y))
     y = y .* y_notnan
 
     u_orig = copy(u)
@@ -43,10 +43,12 @@ function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
     Id_x1                 = spdiagm(x1_deterministic)
 
     if any(x_deterministic)
-        OQ0, OQp = I_nx[find(x_deterministic), :], I_nx[find(!x_deterministic), :]
-        Id(M_t::Array{Int,2}) = spdiagm(OQ0' * all((OQ0 * M_t * OQp') .== 0, 2)[:])
+        OQ0, OQp = I_nx[find(x_deterministic), :], I_nx[find((!).(x_deterministic)), :]
+        # Id(M_t::Array{Int,2}) = spdiagm(OQ0' * all((OQ0 * M_t * OQp') .== 0, 2)[:])
+        Id = M_t -> spdiagm(OQ0' * all((OQ0 * M_t * OQp') .== 0, 2)[:])
     else
-        Id(_::Array{Int,2}) = I0_nx
+        # Id(_::Array{Int,2}) = I0_nx
+        Id = t -> I0_nx
     end
 
     phi(t)  = (pmodel.G(t)' * pmodel.G(t)) \ pmodel.G(t)'
@@ -85,7 +87,7 @@ function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
         =#
 
         if estimate_A || estimate_B || estimate_Q || estimate_x1
-            Qinv    = inv(pmodel.Q(params.Q))  
+            Qinv    = inv(pmodel.Q(params.Q))
             Vinv_   = all_x_deterministic ? I0_nx : phi(1)' * Qinv * phi(1)
         end #if ABQ
 
@@ -104,7 +106,7 @@ function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
 
         function N(t::Int)
             HRHt = HRH(t)
-            O = I_ny[find(y_notnan[:, t] & HRH_nonzero_rows(t)), :]
+            O = I_ny[find(y_notnan[:, t] .& HRH_nonzero_rows(t)), :]
             return I - HRHt * O' * inv(O * HRHt * O') * O
         end #N
 
@@ -221,7 +223,7 @@ function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
                     end #if C
 
                     if estimate_R
-                        I2 = diagm(!y_notnan[:, t])
+                        I2 = diagm(.!y_notnan[:, t])
                         eyy = I2 * (Nt * HRH(t) + Nt * Ct * Pt * Ct' * Nt') * I2 + ey * ey'
                         R_S1 += pmodel.R.D' * pmodel.R.D
                         R_S2 += pmodel.R.D' * vec(xi(t) * (eyy - eyx * Ct' -
@@ -324,7 +326,7 @@ function fit{T}(y::Array{T}, pmodel::ParametrizedSSM, params::SSMParameters;
     ll, ll_prev = Inf, Inf
 
     while (niter > 0) && (fraction_change(ll_prev, ll) > eps)
-        ll_prev = ll 
+        ll_prev = ll
         ll      = em_kernel!(params)
         niter -= 1
     end #while