Merge pull request #10 from fipelle/dev

Dev: new trend parametrisation for DFMs

Project.toml

@@ -1,6 +1,6 @@
 name = "MessyTimeSeriesOptim"
 uuid = "7115d47b-e118-4204-bc33-dbd7ed133e66"
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -10,6 +10,8 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 MessyTimeSeries = "2a88db5c-15f1-4b3e-a070-d1159e8d76cc"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
@@ -25,5 +27,7 @@ test = ["Test"]
 Distributions = "0.25"
 LoopVectorization = "0.12"
 MessyTimeSeries = "0.2"
+Optimization = "3.7"
+OptimizationNLopt = "0.1"
 StableRNGs = "1.0"
 julia = "1.6"

src/MessyTimeSeriesOptim.jl

@@ -5,8 +5,8 @@ module MessyTimeSeriesOptim
 	# Libraries
 	using Dates, Distributed, Logging, LoopVectorization;
 	using Distributions, LinearAlgebra, SparseArrays, StableRNGs, Statistics;
-	using MessyTimeSeries;
-
+	using MessyTimeSeries, Optimization, OptimizationNLopt;
+	
 	# Aliases for MessyTimeSeries
 	find_observed_data = MessyTimeSeries.find_observed_data;
 	update_smoothing_factors! = MessyTimeSeries.update_smoothing_factors!;
@@ -15,15 +15,17 @@ module MessyTimeSeriesOptim
 	# Custom dependencies
 	local_path = dirname(@__FILE__);
 	include("$(local_path)/types.jl");
+	include("$(local_path)/initialisation.jl");
+	include("$(local_path)/estimation.jl");
+	include("$(local_path)/validation.jl");
 	include("$(local_path)/models/dfm.jl");
 	include("$(local_path)/models/var.jl");
 	include("$(local_path)/models/vma.jl");
-	include("$(local_path)/estimation.jl");
-	include("$(local_path)/validation.jl");
 
 	# Export
 	export EstimSettings, DFMSettings, VARSettings, VMASettings, ValidationSettings, HyperGrid;
 	export build_Γ;
+	export initial_univariate_decomposition_kitagawa, initial_univariate_decomposition_llt;
 	export ecm;
 	export select_hyperparameters, fc_err, jackknife_err;
 end

src/estimation.jl

@@ -243,11 +243,12 @@ Initialise M, N and O.
 Initialise M, N and O to nothing.
 """
 function initialise_ecm_stats_measurement(estim::EstimSettings, coordinates_measurement_states::IntVector)
-    M = zeros(estim.n, estim.m);
+    len_measurement_states = length(coordinates_measurement_states);
+    M = zeros(estim.n, len_measurement_states);
     N = Vector{SparseMatrixCSC{Float64, Int64}}(undef, estim.T);
     O = Array{VectorsArray{Float64},1}(undef, estim.T);
     buffer_M = zeros(estim.n);
-    buffer_O = Symmetric(zeros(estim.m, estim.m));
+    buffer_O = Symmetric(zeros(len_measurement_states, len_measurement_states));
     return M, N, O, buffer_M, buffer_O;
 end
 
@@ -294,7 +295,7 @@ Do not update ECM statistics for the measurement equation of models for which th
 """
 function update_ecm_stats_measurement!(barrier_M::FloatMatrix, estim::EstimSettings, smoother_arrays::SmootherArrays, coordinates_measurement_states::IntVector, ind_not_missings::IntVector, t::Int64, Xs::FloatVector, Ps::SymMatrix)
 
-    # Views
+    # Views used for computing M and O (based on `coordinates_measurement_states`)
     Y_obs = @view estim.Y[ind_not_missings, t];
     Xs_view = @view Xs[coordinates_measurement_states];
     Ps_view = @view Ps[coordinates_measurement_states, coordinates_measurement_states];
@@ -310,7 +311,7 @@ function update_ecm_stats_measurement!(barrier_M::FloatMatrix, estim::EstimSetti
     mul!(smoother_arrays.M, smoother_arrays.buffer_M, Xs_view', 1.0, 1.0);
 
     # Update ECM statistics: compute and store N_t
-    smoother_arrays.N[t] = A_transpose*A_transpose';
+    smoother_arrays.N[t] = A_transpose*A_transpose'; # TBD: change the generation of N_{t} -> this is a diagonal matrix with ones at the entries corresponding to series with observed values at time t
 
     # Update ECM statistics: compute O_t
     mul!(smoother_arrays.buffer_O.data, Xs_view, Xs_view');
@@ -480,9 +481,11 @@ function cm_step_X0_P0!(sspace::KalmanSettings, status::SizedKalmanStatus, smoot
     # CM-step for sspace.P0
     mul!(status.online_status.buffer_J2, sspace.P0, Array(smoother_arrays.J2));         # Array(...) helps speeding up mul!(...) while keeping J2 symmetric (i.e., good trade-off wrt other options incl. ".data")
     mul!(status.online_status.buffer_m_m, status.online_status.buffer_J2, sspace.P0);
+    P0_raw = Symmetric(sspace.P0.data - status.online_status.buffer_m_m);
     for ij in coordinates_transition_P0
-        sspace.P0.data[ij] -= status.online_status.buffer_m_m[ij];
+        sspace.P0.data[ij] = P0_raw[ij]; # TBD: symmetry could be exploited in the indexing
     end
+    return P0_raw;
 end
 
 """
@@ -533,10 +536,10 @@ function ksmoother_ecm!(estim::EstimSettings, sspace::KalmanSettings, status::Si
     call_update_smoothing_factors!(sspace, status, smoother_arrays, nothing);
 
     # CM-step for X0 and P0
-    cm_step_X0_P0!(sspace, status, smoother_arrays, coordinates_transition_P0);
+    P0_raw = cm_step_X0_P0!(sspace, status, smoother_arrays, coordinates_transition_P0);
 
     # Finalise calculations of the ECM statistics for the transition equation
-    call_update_ecm_stats_transition!(estim, smoother_arrays, sspace.X0, sspace.P0, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
+    call_update_ecm_stats_transition!(estim, smoother_arrays, sspace.X0, P0_raw, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
 end
 
 #=