From e70232494d8d7987ae4c8064c746e95fa22dbb67 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sat, 23 Dec 2023 20:27:29 +0000 Subject: [PATCH 01/14] Further simplifications in init --- src/initialisation.jl | 147 ++++++++++++++++++++++-------------------- 1 file changed, 78 insertions(+), 69 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 5a2c83a..00ad398 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -1,54 +1,16 @@ """ - update_sspace_B_from_params!(constrained_params::AbstractVector{Float64}, coordinates_free_params_B::BitMatrix, sspace::KalmanSettings) - -Update free coordinates in `sspace.B` from `constrained_params`. -""" -function update_sspace_B_from_params!(constrained_params::AbstractVector{Float64}, coordinates_free_params_B::BitMatrix, sspace::KalmanSettings) - sspace.B[coordinates_free_params_B] .= constrained_params[1:sum(coordinates_free_params_B)]; -end - -""" - update_sspace_Q_from_params!(constrained_params::AbstractVector{Float64}, coordinates_free_params_B::BitMatrix, sspace::KalmanSettings) - -Update `sspace.Q` from `constrained_params`. -""" -function update_sspace_Q_from_params!(constrained_params::AbstractVector{Float64}, coordinates_free_params_B::BitMatrix, sspace::KalmanSettings) - - # Find relevant coordinates - n_series = size(sspace.B, 1); - n_trends_ext = findall(sspace.B[1,:] .== 1)[2]-1; - n_trends = n_trends_ext/2 |> Int64; - - # Trends dictionary - trends_dict = Dict(i=>j for (j, i) in enumerate(1:2:n_trends_ext)); - - # Break down parameters - params_ratios = constrained_params[sum(coordinates_free_params_B)+1:sum(coordinates_free_params_B)+1+n_series]; - params_variances = constrained_params[sum(coordinates_free_params_B)+2+n_series:end]; - - # Ratios to variances - params_ratios_as_variances = copy(params_ratios); - params_ratios_as_variances[end] *= params_variances[1]; # first common cycle - for i=1:length(params_ratios_as_variances)-1 - coordinates_trends_in_B = findall(sspace.B[i, 1:n_trends_ext] .!= 0.0); - coordinates_trends_in_params = getindex.(Ref(trends_dict), coordinates_trends_in_B); - params_ratios_as_variances[i] *= sum(view(params_variances, coordinates_trends_in_params)); # idio cycles - end - - # Merge variances - merged_params_variances = vcat(params_variances[1:n_trends], params_ratios_as_variances, params_variances[n_trends+1:end]); - - # Update sspace.Q - sspace.Q.data[diagind(sspace.Q.data)] .= merged_params_variances; -end - -""" - update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0::BitMatrix, sspace::KalmanSettings) + update_sspace_DQD_and_P0_from_params!( + sspace :: KalmanSettings, + coordinates_free_params_P0 :: BitMatrix + ) Update `sspace.DQD` and the free entries of `sspace.P0`. """ -function update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0::BitMatrix, sspace::KalmanSettings) - +function update_sspace_DQD_and_P0_from_params!( + sspace :: KalmanSettings, + coordinates_free_params_P0 :: BitMatrix +) + # Update `sspace.DQD` sspace.DQD.data .= Symmetric(sspace.D*sspace.Q*sspace.D').data; @@ -65,20 +27,50 @@ function update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0::BitMa end """ - fmin!(constrained_params::AbstractVector{Float64}, sspace::KalmanSettings, coordinates_free_params_B::BitMatrix, coordinates_free_params_P0::BitMatrix) + fmin!( + constrained_params :: AbstractVector{Float64}, + sspace :: KalmanSettings, + coordinates_free_params_B :: BitMatrix, + coordinates_free_params_Q :: BitMatrix, + coordinates_free_params_P0 :: BitMatrix + ) Return -1 times the log-likelihood function (or a large number in case of numerical problems). """ -function fmin!(constrained_params::AbstractVector{Float64}, sspace::KalmanSettings, coordinates_free_params_B::BitMatrix, coordinates_free_params_P0::BitMatrix) +function fmin!( + constrained_params :: AbstractVector{Float64}, + sspace :: KalmanSettings, + coordinates_free_params_B :: BitMatrix, + coordinates_free_params_Q :: BitMatrix, + coordinates_free_params_P0 :: BitMatrix +) + + # Useful shortcut + counter = 0; + + # Free parameters in the measurement equation coefficients + if sum(coordinates_free_params_B) > 0 + sspace.B[coordinates_free_params_B] .= constrained_params[1:sum(coordinates_free_params_B)]; + counter += sum(coordinates_free_params_B); + end - # Update sspace accordingly - update_sspace_B_from_params!(constrained_params, coordinates_free_params_B, sspace); - update_sspace_Q_from_params!(constrained_params, coordinates_free_params_B, sspace); + # Free parameters in the state equation coefficients + #= + if sum(coordinates_free_params_C) > 0 + sspace.C[coordinates_free_params_C] .= constrained_params[counter+1:counter+sum(coordinates_free_params_C)]; + counter += sum(coordinates_free_params_C); + end + =# + + # Free parameters in the state equation covariance matrix + if sum(coordinates_free_params_Q) > 0 + sspace.Q.data[coordinates_free_params_Q] .= constrained_params[counter+1:counter+end]; + end # Determine whether Q is problematic if (sum(isnan.(sspace.Q)) == 0) && (sum(isinf.(sspace.Q)) == 0) is_Q_non_problematic = true; - update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0, sspace); + update_sspace_DQD_and_P0_from_params!(sspace, coordinates_free_params_P0); else is_Q_non_problematic = false; end @@ -121,13 +113,21 @@ APIs to call `fmin!` with Tuple parameters. call_fmin!(constrained_params::AbstractVector{Float64}, tuple_fmin_args::Tuple) = fmin!(constrained_params, tuple_fmin_args...); """ - initial_sspace_structure(data::Union{FloatMatrix, JMatrix{Float64}}, estim::EstimSettings; first_step::Bool=false) + initial_sspace_structure( + data :: Union{FloatMatrix, JMatrix{Float64}}, + estim :: EstimSettings; + first_step :: Bool = false + ) Get initial state-space parameters and relevant coordinates. Trends are modelled using the Kitagawa representation. """ -function initial_sspace_structure(data::Union{FloatMatrix, JMatrix{Float64}}, estim::EstimSettings; first_step::Bool=false) +function initial_sspace_structure( + data :: Union{FloatMatrix, JMatrix{Float64}}, + estim :: EstimSettings; + first_step :: Bool = false +) # `n` for initialisation (it may differ from the one in `estim`) n_series_in_data = size(data, 1); @@ -151,6 +151,9 @@ function initial_sspace_structure(data::Union{FloatMatrix, JMatrix{Float64}}, es # Setup loading matrix B = hcat(B_trends, B_idio_cycles, B_common_cycles); + + # Setup BitMatrix for the free parameters in `B` + # NOTE: shortcut to update the factor loadings (not if first_step) coordinates_free_params_B = B .> 1.0; # Setup covariance matrix measurement error @@ -176,8 +179,13 @@ function initial_sspace_structure(data::Union{FloatMatrix, JMatrix{Float64}}, es # Setup covariance matrix of the states' innovation # NOTE: all diagonal elements are estimated during the initialisation - Q = Symmetric(Matrix(1.0I, estim.n_trends + n_series_in_data + estim.n_cycles, estim.n_trends + n_series_in_data + estim.n_cycles)); + Q = Symmetric(Matrix(1e-4*I, estim.n_trends + n_series_in_data + estim.n_cycles, estim.n_trends + n_series_in_data + estim.n_cycles)); + # Setup BitMatrix for the free parameters in `Q` + # NOTE: All diagonal elements but those referring to the trends, whose variances are fixed to 1e-4 + coordinates_free_params_Q = Q .== 1; + coordinates_free_params_Q[1:estim.n_trends, 1:estim.n_trends] .= false; + # Setup selection matrix D for the trends D_trends_template = [1.0; 0.0]; D_trends = cat(dims=[1,2], [D_trends_template for i in 1:estim.n_trends]...); @@ -225,11 +233,13 @@ function initial_sspace_structure(data::Union{FloatMatrix, JMatrix{Float64}}, es # Setup initial conditions X0 = vcat(X0_trends, X0_idio_cycles, X0_common_cycles); P0 = Symmetric(cat(dims=[1,2], P0_trends, P0_idio_cycles, P0_common_cycles)); + + # Setup BitMatrix for the free parameters in `P0` coordinates_free_params_P0 = P0 .!= 0.0; coordinates_free_params_P0[1:2*estim.n_trends, 1:2*estim.n_trends] .= false; # Return state-space matrices and relevant coordinates - return B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_P0; + return B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0; end """ @@ -240,25 +250,24 @@ Run first step of the initialisation to find reasonable initial guesses. function initial_detrending_step_1(Y_trimmed::JMatrix{Float64}, estim::EstimSettings, n_trimmed::Int64) # Get initial state-space parameters and relevant coordinates - B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim, first_step=true); + B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim, first_step=true); # Set KalmanSettings sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); # Initial guess for the parameters - params_0 = vcat( - 1e+4*ones(1+n_trimmed), - 1e-3*ones(estim.n_trends), - ); - params_lb = vcat(1e+2*ones(1+n_trimmed), 1e-6*ones(estim.n_trends)); - params_ub = vcat(1e+6*ones(1+n_trimmed), ones(estim.n_trends)); + params_0 = 1e-1*ones(1+n_trimmed); # variances of the cycles' innovations + params_lb = 1e-3*ones(1+n_trimmed); # NOTE: > 1e-4 + params_ub = 1e+3*ones(1+n_trimmed); # Maximum likelihood - tuple_fmin_args = (sspace, coordinates_free_params_B, coordinates_free_params_P0); + tuple_fmin_args = (sspace, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0); prob = OptimizationFunction(call_fmin!) prob = OptimizationProblem(prob, params_0, tuple_fmin_args, lb=params_lb, ub=params_ub); - res_optim = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2); - + + # Recover optimum + res_optim = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2).u; + # Update `sspace` from `res_optim` update_sspace_Q_from_params!(res_optim.u, coordinates_free_params_B, sspace); update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0, sspace); @@ -392,7 +401,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e params_0, _, smoothed_cycles_0 = MessyTimeSeriesOptim.initial_detrending_step_1(Y_trimmed, estim, n_trimmed); # Get initial state-space parameters and relevant coordinates - B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_P0 = MessyTimeSeriesOptim.initial_sspace_structure(Y_trimmed, estim); + B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = MessyTimeSeriesOptim.initial_sspace_structure(Y_trimmed, estim); # Set `coordinates_free_params_B` to `false` coordinates_free_params_B .= false; @@ -462,7 +471,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Maximum likelihood @info("Initialisation > running optimizer") - tuple_fmin_args = (sspace, coordinates_free_params_B, coordinates_free_params_P0); + tuple_fmin_args = (sspace, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0); prob = OptimizationFunction(call_fmin!) prob = OptimizationProblem(prob, params_0, tuple_fmin_args, lb=params_lb, ub=params_ub); res_optim = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2); From 2f62da0d7654088a89a0d6cd7a6a44a9252becf1 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sat, 23 Dec 2023 20:32:15 +0000 Subject: [PATCH 02/14] Added shortcut for sspace update --- src/initialisation.jl | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 00ad398..9cbcff9 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -27,7 +27,7 @@ function update_sspace_DQD_and_P0_from_params!( end """ - fmin!( + update_sspace_from_params!( constrained_params :: AbstractVector{Float64}, sspace :: KalmanSettings, coordinates_free_params_B :: BitMatrix, @@ -35,9 +35,9 @@ end coordinates_free_params_P0 :: BitMatrix ) -Return -1 times the log-likelihood function (or a large number in case of numerical problems). +Update free parameters in `sspace`. """ -function fmin!( +function update_sspace_from_params!( constrained_params :: AbstractVector{Float64}, sspace :: KalmanSettings, coordinates_free_params_B :: BitMatrix, @@ -81,6 +81,37 @@ function fmin!( else is_P0_non_problematic = false; end + + return is_Q_non_problematic, is_P0_non_problematic; +end + +""" + fmin!( + constrained_params :: AbstractVector{Float64}, + sspace :: KalmanSettings, + coordinates_free_params_B :: BitMatrix, + coordinates_free_params_Q :: BitMatrix, + coordinates_free_params_P0 :: BitMatrix + ) + +Return -1 times the log-likelihood function (or a large number in case of numerical problems). +""" +function fmin!( + constrained_params :: AbstractVector{Float64}, + sspace :: KalmanSettings, + coordinates_free_params_B :: BitMatrix, + coordinates_free_params_Q :: BitMatrix, + coordinates_free_params_P0 :: BitMatrix +) + + # Update `sspace` free parameters + is_Q_non_problematic, is_P0_non_problematic = update_sspace_from_params!( + constrained_params, + sspace, + coordinates_free_params_B, + coordinates_free_params_Q, + coordinates_free_params_P0 + ) # Regular run if is_Q_non_problematic && is_P0_non_problematic From baccc374f7adb95342c5b02cb19fd365399a8317 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sat, 23 Dec 2023 20:37:31 +0000 Subject: [PATCH 03/14] Ultimated simplification --- src/initialisation.jl | 47 ++++++++++++++++++++++++++----------------- 1 file changed, 29 insertions(+), 18 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 9cbcff9..bab1425 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -111,7 +111,7 @@ function fmin!( coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 - ) + ); # Regular run if is_Q_non_problematic && is_P0_non_problematic @@ -297,11 +297,16 @@ function initial_detrending_step_1(Y_trimmed::JMatrix{Float64}, estim::EstimSett prob = OptimizationProblem(prob, params_0, tuple_fmin_args, lb=params_lb, ub=params_ub); # Recover optimum - res_optim = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2).u; + constrained_params = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2).u; - # Update `sspace` from `res_optim` - update_sspace_Q_from_params!(res_optim.u, coordinates_free_params_B, sspace); - update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0, sspace); + # Update `sspace` free parameters + _ = update_sspace_from_params!( + constrained_params, + sspace, + coordinates_free_params_B, + coordinates_free_params_Q, + coordinates_free_params_P0 + ); # Recover smoothed states status = kfilter_full_sample(sspace); @@ -318,7 +323,7 @@ function initial_detrending_step_1(Y_trimmed::JMatrix{Float64}, estim::EstimSett end # Return minimizer - return res_optim.u, smoothed_states, smoothed_cycles; + return constrained_params, smoothed_states, smoothed_cycles; end """ @@ -429,10 +434,10 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Recover initial guess from step 1 @info("Initialisation > warm start") - params_0, _, smoothed_cycles_0 = MessyTimeSeriesOptim.initial_detrending_step_1(Y_trimmed, estim, n_trimmed); + _, _, smoothed_cycles_0 = MessyTimeSeriesOptim.initial_detrending_step_1(Y_trimmed, estim, n_trimmed); # Get initial state-space parameters and relevant coordinates - B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = MessyTimeSeriesOptim.initial_sspace_structure(Y_trimmed, estim); + B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim); # Set `coordinates_free_params_B` to `false` coordinates_free_params_B .= false; @@ -495,23 +500,29 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Set KalmanSettings sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); - # Update `params_0` - # The variance for the cycles is re-calibrated more aggressively given that in the warm start there are no common cycles - params_lb = vcat(1e+2*ones(1+n_trimmed), params_0[2+n_trimmed:end]/10); - params_ub = vcat(1e+6*ones(1+n_trimmed), params_0[2+n_trimmed:end]*10); - + # Initial guess for the parameters + params_0 = 1e-1*ones(1+n_trimmed); # variances of the cycles' innovations + params_lb = 1e-3*ones(1+n_trimmed); # NOTE: > 1e-4 + params_ub = 1e+3*ones(1+n_trimmed); + # Maximum likelihood @info("Initialisation > running optimizer") tuple_fmin_args = (sspace, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0); prob = OptimizationFunction(call_fmin!) prob = OptimizationProblem(prob, params_0, tuple_fmin_args, lb=params_lb, ub=params_ub); - res_optim = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2); - # Update sspace accordingly - update_sspace_B_from_params!(res_optim.u, coordinates_free_params_B, sspace); - update_sspace_Q_from_params!(res_optim.u, coordinates_free_params_B, sspace); - update_sspace_DQD_and_P0_from_params!(coordinates_free_params_P0, sspace); + # Recover optimum + constrained_params = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2).u; + # Update `sspace` free parameters + _ = update_sspace_from_params!( + constrained_params, + sspace, + coordinates_free_params_B, + coordinates_free_params_Q, + coordinates_free_params_P0 + ); + # Recover smoothed states status = kfilter_full_sample(sspace); smoothed_states_container, _ = ksmoother(sspace, status); From 972b8bd19e2f0144835a212e0d24ca1025920990 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sun, 24 Dec 2023 02:28:57 +0000 Subject: [PATCH 04/14] minor fix --- src/initialisation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index bab1425..09931b6 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -64,7 +64,7 @@ function update_sspace_from_params!( # Free parameters in the state equation covariance matrix if sum(coordinates_free_params_Q) > 0 - sspace.Q.data[coordinates_free_params_Q] .= constrained_params[counter+1:counter+end]; + sspace.Q.data[coordinates_free_params_Q] .= constrained_params[counter+1:end]; end # Determine whether Q is problematic @@ -214,7 +214,7 @@ function initial_sspace_structure( # Setup BitMatrix for the free parameters in `Q` # NOTE: All diagonal elements but those referring to the trends, whose variances are fixed to 1e-4 - coordinates_free_params_Q = Q .== 1; + coordinates_free_params_Q = Q .!= 0.0; coordinates_free_params_Q[1:estim.n_trends, 1:estim.n_trends] .= false; # Setup selection matrix D for the trends From c0806a240f30190ac2d4876c9676307bbd5b5b0a Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sun, 24 Dec 2023 17:55:52 +0000 Subject: [PATCH 05/14] Deprecated first step --- src/initialisation.jl | 22 ++++++---------------- 1 file changed, 6 insertions(+), 16 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 09931b6..e8f41b8 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -145,9 +145,8 @@ call_fmin!(constrained_params::AbstractVector{Float64}, tuple_fmin_args::Tuple) """ initial_sspace_structure( - data :: Union{FloatMatrix, JMatrix{Float64}}, - estim :: EstimSettings; - first_step :: Bool = false + data :: Union{FloatMatrix, JMatrix{Float64}}, + estim :: EstimSettings ) Get initial state-space parameters and relevant coordinates. @@ -155,9 +154,8 @@ Get initial state-space parameters and relevant coordinates. Trends are modelled using the Kitagawa representation. """ function initial_sspace_structure( - data :: Union{FloatMatrix, JMatrix{Float64}}, - estim :: EstimSettings; - first_step :: Bool = false + data :: Union{FloatMatrix, JMatrix{Float64}}, + estim :: EstimSettings ) # `n` for initialisation (it may differ from the one in `estim`) @@ -176,15 +174,11 @@ function initial_sspace_structure( B_common_cycles[1, i] = 1.0; end - if first_step - B_common_cycles[B_common_cycles .!= 1.0] .= 0.0; # no common cycles - end - # Setup loading matrix B = hcat(B_trends, B_idio_cycles, B_common_cycles); # Setup BitMatrix for the free parameters in `B` - # NOTE: shortcut to update the factor loadings (not if first_step) + # NOTE: shortcut to update the factor loadings coordinates_free_params_B = B .> 1.0; # Setup covariance matrix measurement error @@ -196,11 +190,7 @@ function initial_sspace_structure( # Setup transition matrix for the idiosyncratic cycles C_idio_cycles = Matrix(0.1I, n_series_in_data, n_series_in_data); - if first_step - C_idio_cycles[1, 1] = 0.0; # set to noise - C_idio_cycles[C_idio_cycles .== 0.1] .*= 9.0; # set to persistent cycles (as persistent as the common cycles) - end - + # Setup transition matrix for the common cycles C_common_cycles_template = companion_form([0.9 zeros(1, estim.lags-1)], extended=false); C_common_cycles = cat(dims=[1,2], [C_common_cycles_template for i in 1:estim.n_cycles]...); From 684344f5efc5fbe2cccf16b01709755941b368d0 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sun, 24 Dec 2023 17:56:25 +0000 Subject: [PATCH 06/14] Deprecated functions for warm start --- src/initialisation.jl | 134 +----------------------------------------- 1 file changed, 1 insertion(+), 133 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index e8f41b8..30f1d06 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -190,7 +190,7 @@ function initial_sspace_structure( # Setup transition matrix for the idiosyncratic cycles C_idio_cycles = Matrix(0.1I, n_series_in_data, n_series_in_data); - + # Setup transition matrix for the common cycles C_common_cycles_template = companion_form([0.9 zeros(1, estim.lags-1)], extended=false); C_common_cycles = cat(dims=[1,2], [C_common_cycles_template for i in 1:estim.n_cycles]...); @@ -263,138 +263,6 @@ function initial_sspace_structure( return B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0; end -""" - initial_detrending_step_1(Y_trimmed::JMatrix{Float64}, estim::EstimSettings, n_trimmed::Int64) - -Run first step of the initialisation to find reasonable initial guesses. -""" -function initial_detrending_step_1(Y_trimmed::JMatrix{Float64}, estim::EstimSettings, n_trimmed::Int64) - - # Get initial state-space parameters and relevant coordinates - B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim, first_step=true); - - # Set KalmanSettings - sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); - - # Initial guess for the parameters - params_0 = 1e-1*ones(1+n_trimmed); # variances of the cycles' innovations - params_lb = 1e-3*ones(1+n_trimmed); # NOTE: > 1e-4 - params_ub = 1e+3*ones(1+n_trimmed); - - # Maximum likelihood - tuple_fmin_args = (sspace, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0); - prob = OptimizationFunction(call_fmin!) - prob = OptimizationProblem(prob, params_0, tuple_fmin_args, lb=params_lb, ub=params_ub); - - # Recover optimum - constrained_params = solve(prob, NLopt.LN_SBPLX, abstol=1e-3, reltol=1e-2).u; - - # Update `sspace` free parameters - _ = update_sspace_from_params!( - constrained_params, - sspace, - coordinates_free_params_B, - coordinates_free_params_Q, - coordinates_free_params_P0 - ); - - # Recover smoothed states - status = kfilter_full_sample(sspace); - smoothed_states_container, _ = ksmoother(sspace, status); - smoothed_states = hcat(smoothed_states_container...); - - # Recover smoothed cycles - smoothed_cycles = smoothed_states[2*estim.n_trends+1:2*estim.n_trends+n_trimmed, :]; - for i=1:n_trimmed - last_state_for_ith_series = findlast(B[i, :] .== 1.0); - if last_state_for_ith_series > 2*estim.n_trends+n_trimmed - smoothed_cycles[i, :] .+= smoothed_states[last_state_for_ith_series, :]; - end - end - - # Return minimizer - return constrained_params, smoothed_states, smoothed_cycles; -end - -""" - initialise_common_cycle(estim::EstimSettings, residual_data::FloatMatrix, coordinates_current_block::IntVector) - -Initialise current common cycle via PCA. -""" -function initialise_common_cycle(estim::EstimSettings, residual_data::FloatMatrix, coordinates_current_block::IntVector) - - # Convenient shortcuts - data_current_block = residual_data[coordinates_current_block, :]; - data_current_block_standardised = standardise(data_current_block); - - # Compute PCA loadings - eigen_val, eigen_vect = eigen(Symmetric(cov(data_current_block_standardised, dims=2))); - loadings = eigen_vect[:, sortperm(-abs.(eigen_val))[1]]; - - # Compute PCA factor - pc1 = permutedims(loadings)*data_current_block_standardised; - - # Rescale PCA loadings to match the original scale - loadings .*= std(data_current_block, dims=2)[:]; - - # Rescale PC1 wrt the first series - pc1 .*= loadings[1]; - loadings ./= loadings[1]; - - if estim.lags > 1 - - # Backcast `pc1` first - - # Reverse `pc1` time order to predict the past (i.e., backcast) - pc1_reversed = reverse(pc1); - - # Estimate ridge backcast coefficients - pc1_reversed_y, pc1_reversed_x = lag(pc1_reversed, estim.lags); - ar_coeff_backcast = pc1_reversed_y*pc1_reversed_x'/Symmetric(pc1_reversed_x*pc1_reversed_x' + estim.Γ); - enforce_causality_and_invertibility!(ar_coeff_backcast); - - # Generate backcast for `pc1` - for t=1:estim.lags - backcast_x = pc1[1:estim.lags]; - pc1 = hcat(ar_coeff_backcast*backcast_x, pc1); - end - - # Lag principal component - pc1_y, pc1_x = lag(pc1, estim.lags); - - # Initialise complete loadings - complete_loadings = zeros(length(loadings), estim.lags); - - # Regress one variable at the time on `pc1` - pc1_x_shifted_with_backcast = vcat(pc1_y, pc1_x[1:end-1, :]); - pc1_x_shifted = pc1_x_shifted_with_backcast[:, estim.lags+1:end]; - for i in axes(data_current_block, 1) - if i == 1 - complete_loadings[1, 1] = 1.0; # identification - else - data_current_block_yi, _ = lag(permutedims(data_current_block[i, :]), estim.lags); - complete_loadings[i, :] = data_current_block_yi*pc1_x_shifted'/Symmetric(pc1_x_shifted*pc1_x_shifted' + estim.Γ); - end - end - - # Estimate autoregressive dynamics - ar_coefficients = pc1_y*pc1_x'/Symmetric(pc1_x*pc1_x' + estim.Γ); - - # Estimate var-cov matrix of the residuals - ar_residuals = pc1_y - ar_coefficients*pc1_x; - ar_residuals_variance = (ar_residuals*ar_residuals')/length(ar_residuals); - - # Explained data - explained_data = complete_loadings*pc1_x_shifted_with_backcast; - - else - error("`estim.lags` must be greater than one!"); - end - - # Return output - return complete_loadings, ar_coefficients, ar_residuals_variance, explained_data; -end - """ initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, estim::EstimSettings) From e9ff50030a8ff50e9c9e8dc0d13120bd319af56a Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sun, 24 Dec 2023 17:59:10 +0000 Subject: [PATCH 07/14] Deprecated warm start --- src/initialisation.jl | 61 +------------------------------------------ 1 file changed, 1 insertion(+), 60 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 30f1d06..10be975 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -289,72 +289,13 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e last_ind = findlast(sum(ismissing.(Y_untrimmed), dims=1) .== 0)[2]; Y_trimmed = Y_untrimmed[:, first_ind:last_ind] |> JMatrix{Float64}; n_trimmed = size(Y_trimmed, 1); - - # Recover initial guess from step 1 - @info("Initialisation > warm start") - _, _, smoothed_cycles_0 = MessyTimeSeriesOptim.initial_detrending_step_1(Y_trimmed, estim, n_trimmed); - + # Get initial state-space parameters and relevant coordinates B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim); # Set `coordinates_free_params_B` to `false` coordinates_free_params_B .= false; - # Determine which series can load on each cycle - boolean_coordinates_blocks = (estim.cycles_skeleton .!= 0) .| estim.cycles_free_params; - coordinates_blocks = [findall(boolean_coordinates_blocks[:, i]) for i=1:estim.n_cycles]; - - # Initialise common cycles iteratively via PCA - residual_data = copy(smoothed_cycles_0); - - # Loop over common cycles - for i=1:estim.n_cycles - - # Pointers - coordinates_current_block = coordinates_blocks[i]; - - # Loadings for the i-th common cycle - complete_loadings, ar_coefficients, ar_residuals_variance, explained_data = initialise_common_cycle(estim, residual_data, coordinates_current_block); - - # Position of the i-th common cycle in the state-space representation - coordinates_pc1 = 1 + (i-1)*estim.lags + 2*estim.n_trends + n_trimmed; - - # Position `complete_loadings` in `B` - B[:, coordinates_pc1:coordinates_pc1+estim.lags-1] = complete_loadings; - - # Position `ar_dynamics` in `C` - C[coordinates_pc1, coordinates_pc1:coordinates_pc1+estim.lags-1] .= ar_coefficients[:]; - - # Position `resid_variance` in `Q` - Q[estim.n_trends + n_trimmed + i, estim.n_trends + n_trimmed + i] = ar_residuals_variance[1]; - - # Recompute residual data - residual_data[coordinates_current_block, :] .-= explained_data; - end - - idio_ar_coefficients = zeros(n_trimmed); - idio_residuals_variance = zeros(n_trimmed); - - for i=1:n_trimmed - - # Lag residual data - idio_y, idio_x = lag(permutedims(residual_data[i, :]), 1); - - # Idiosyncratic cycles coefficients - idio_ar_coefficient = idio_y*idio_x'/Symmetric(idio_x*idio_x'); - idio_ar_coefficients[i] = idio_ar_coefficient[1]; - - # Idiosyncratic cycles variance - final_residuals = idio_y - idio_ar_coefficient*idio_x; - final_residuals_variance = (final_residuals*final_residuals')/length(final_residuals); - idio_residuals_variance[i] = final_residuals_variance[1]; - end - - C_idio_view = @view C[2*estim.n_trends+1:2*estim.n_trends+n_trimmed, 2*estim.n_trends+1:2*estim.n_trends+n_trimmed]; - C_idio_view[diagind(C_idio_view)] .= idio_ar_coefficients; - Q_idio_view = @view Q[estim.n_trends+1:estim.n_trends+n_trimmed, estim.n_trends+1:estim.n_trends+n_trimmed]; - Q_idio_view[diagind(Q_idio_view)] .= idio_residuals_variance; - # Set KalmanSettings sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); From a4c9b982ef450953eefd893aada3c86d3c8ec7e2 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sun, 24 Dec 2023 18:03:27 +0000 Subject: [PATCH 08/14] Fixed initial guess --- src/initialisation.jl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 10be975..33187fc 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -289,7 +289,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e last_ind = findlast(sum(ismissing.(Y_untrimmed), dims=1) .== 0)[2]; Y_trimmed = Y_untrimmed[:, first_ind:last_ind] |> JMatrix{Float64}; n_trimmed = size(Y_trimmed, 1); - + # Get initial state-space parameters and relevant coordinates B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim); @@ -300,9 +300,22 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); # Initial guess for the parameters - params_0 = 1e-1*ones(1+n_trimmed); # variances of the cycles' innovations - params_lb = 1e-3*ones(1+n_trimmed); # NOTE: > 1e-4 - params_ub = 1e+3*ones(1+n_trimmed); + params_0 = vcat( + ones(sum(coordinates_free_params_B)), # loadings + 1e-1*ones(1+n_trimmed); # variances of the cycles' innovations + ) + + # Lower bound + params_lb = vcat( + -10*ones(sum(coordinates_free_params_B)), # loadings + 1e-2*ones(1+n_trimmed); # variances of the cycles' innovations (NOTE: > 1e-4) + ); + + # Upper bound + params_ub = vcat( + +10*ones(sum(coordinates_free_params_B)), # loadings + 1e+2*ones(1+n_trimmed); # variances of the cycles' innovations + ); # Maximum likelihood @info("Initialisation > running optimizer") From bb05d173cfbdc7e736a6da17f4fd426985266cb3 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Mon, 25 Dec 2023 17:36:01 +0000 Subject: [PATCH 09/14] estimate trends var too --- src/initialisation.jl | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 33187fc..569b540 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -64,7 +64,23 @@ function update_sspace_from_params!( # Free parameters in the state equation covariance matrix if sum(coordinates_free_params_Q) > 0 - sspace.Q.data[coordinates_free_params_Q] .= constrained_params[counter+1:end]; + + # Update cycles' innovations variances + sspace.Q.data[coordinates_free_params_Q] .= constrained_params[counter+1:counter+sum(coordinates_free_params_Q)]; + + # Update counter + counter += sum(coordinates_free_params_Q); + + # Trends' innovations variance ratios + var_ratios = @view constrained_params[counter+1:end]; + + # Update trends' innovations variances + n_trends = length(var_ratios); + for i=1:n_trends + coordinates_series_on_ith_trend = sspace.B[:, 1+(i-1)*2] .== 1; + var_cycles_innovations = sspace.B[coordinates_series_on_ith_trend, :]*sspace.D*sspace.Q*[zeros(n_trends); ones(size(sspace.Q, 1)-n_trends)]; + sspace.Q.data[i, i] = minimum(var_cycles_innovations)*var_ratios[i]; + end end # Determine whether Q is problematic @@ -190,7 +206,8 @@ function initial_sspace_structure( # Setup transition matrix for the idiosyncratic cycles C_idio_cycles = Matrix(0.1I, n_series_in_data, n_series_in_data); - + C_idio_cycles[1, 1] = 0.0; + # Setup transition matrix for the common cycles C_common_cycles_template = companion_form([0.9 zeros(1, estim.lags-1)], extended=false); C_common_cycles = cat(dims=[1,2], [C_common_cycles_template for i in 1:estim.n_cycles]...); @@ -294,7 +311,8 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim); # Set `coordinates_free_params_B` to `false` - coordinates_free_params_B .= false; + coordinates_free_params_B[:, end-2:end] .= false; + B[:, end-2:end] .= 0.0; # Set KalmanSettings sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); @@ -302,19 +320,22 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Initial guess for the parameters params_0 = vcat( ones(sum(coordinates_free_params_B)), # loadings - 1e-1*ones(1+n_trimmed); # variances of the cycles' innovations + 1e-1*ones(1+n_trimmed), # variances of the cycles' innovations + 1e-3*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ) - + # Lower bound params_lb = vcat( -10*ones(sum(coordinates_free_params_B)), # loadings - 1e-2*ones(1+n_trimmed); # variances of the cycles' innovations (NOTE: > 1e-4) + 1e-3*ones(1+n_trimmed), # variances of the cycles' innovations + 1e-6*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); # Upper bound params_ub = vcat( +10*ones(sum(coordinates_free_params_B)), # loadings - 1e+2*ones(1+n_trimmed); # variances of the cycles' innovations + 1e+3*ones(1+n_trimmed), # variances of the cycles' innovations + ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); # Maximum likelihood From afb85ead808e63e3f1f31eacbf282a443491c826 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Tue, 26 Dec 2023 01:07:02 +0000 Subject: [PATCH 10/14] Aligned with older version --- src/initialisation.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 569b540..906cafa 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -78,7 +78,7 @@ function update_sspace_from_params!( n_trends = length(var_ratios); for i=1:n_trends coordinates_series_on_ith_trend = sspace.B[:, 1+(i-1)*2] .== 1; - var_cycles_innovations = sspace.B[coordinates_series_on_ith_trend, :]*sspace.D*sspace.Q*[zeros(n_trends); ones(size(sspace.Q, 1)-n_trends)]; + var_cycles_innovations = (sspace.B[coordinates_series_on_ith_trend, :] .!= 0.0)*sspace.D*sspace.Q*[zeros(n_trends); ones(size(sspace.Q, 1)-n_trends)]; sspace.Q.data[i, i] = minimum(var_cycles_innovations)*var_ratios[i]; end end @@ -320,22 +320,22 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Initial guess for the parameters params_0 = vcat( ones(sum(coordinates_free_params_B)), # loadings - 1e-1*ones(1+n_trimmed), # variances of the cycles' innovations - 1e-3*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) + 1e+1*ones(1+n_trimmed), # variances of the cycles' innovations + 1e-4*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ) # Lower bound params_lb = vcat( -10*ones(sum(coordinates_free_params_B)), # loadings - 1e-3*ones(1+n_trimmed), # variances of the cycles' innovations + 1e-4*ones(1+n_trimmed), # variances of the cycles' innovations 1e-6*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); # Upper bound params_ub = vcat( +10*ones(sum(coordinates_free_params_B)), # loadings - 1e+3*ones(1+n_trimmed), # variances of the cycles' innovations - ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) + 1e+6*ones(1+n_trimmed), # variances of the cycles' innovations + 1e-2*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); # Maximum likelihood From 0e6d97893cdeb8d4a799114ba918f2105f80ee7d Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Tue, 26 Dec 2023 02:29:33 +0000 Subject: [PATCH 11/14] Less extreme bounds --- src/initialisation.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 906cafa..64a0ace 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -207,7 +207,7 @@ function initial_sspace_structure( # Setup transition matrix for the idiosyncratic cycles C_idio_cycles = Matrix(0.1I, n_series_in_data, n_series_in_data); C_idio_cycles[1, 1] = 0.0; - + # Setup transition matrix for the common cycles C_common_cycles_template = companion_form([0.9 zeros(1, estim.lags-1)], extended=false); C_common_cycles = cat(dims=[1,2], [C_common_cycles_template for i in 1:estim.n_cycles]...); @@ -320,7 +320,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Initial guess for the parameters params_0 = vcat( ones(sum(coordinates_free_params_B)), # loadings - 1e+1*ones(1+n_trimmed), # variances of the cycles' innovations + ones(1+n_trimmed), # variances of the cycles' innovations 1e-4*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ) @@ -334,7 +334,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Upper bound params_ub = vcat( +10*ones(sum(coordinates_free_params_B)), # loadings - 1e+6*ones(1+n_trimmed), # variances of the cycles' innovations + 1e+4*ones(1+n_trimmed), # variances of the cycles' innovations 1e-2*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); From 6ff60b41f851be7539a48955f01a76311bd8d127 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sat, 20 Jan 2024 16:39:33 +0000 Subject: [PATCH 12/14] Finetuning for variance contraints --- src/initialisation.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 64a0ace..0684643 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -78,7 +78,7 @@ function update_sspace_from_params!( n_trends = length(var_ratios); for i=1:n_trends coordinates_series_on_ith_trend = sspace.B[:, 1+(i-1)*2] .== 1; - var_cycles_innovations = (sspace.B[coordinates_series_on_ith_trend, :] .!= 0.0)*sspace.D*sspace.Q*[zeros(n_trends); ones(size(sspace.Q, 1)-n_trends)]; + var_cycles_innovations = sspace.B[coordinates_series_on_ith_trend, :]*sspace.D*sspace.Q*[zeros(n_trends); ones(size(sspace.Q, 1)-n_trends)]; sspace.Q.data[i, i] = minimum(var_cycles_innovations)*var_ratios[i]; end end @@ -311,8 +311,8 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e B, R, C, D, Q, X0, P0, coordinates_free_params_B, coordinates_free_params_Q, coordinates_free_params_P0 = initial_sspace_structure(Y_trimmed, estim); # Set `coordinates_free_params_B` to `false` - coordinates_free_params_B[:, end-2:end] .= false; - B[:, end-2:end] .= 0.0; + #coordinates_free_params_B[:, end-2:end] .= false; + #B[:, end-2:end] .= 0.0; # Set KalmanSettings sspace = KalmanSettings(Y_trimmed, B, R, C, D, Q, X0, P0, compute_loglik=true); @@ -321,7 +321,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e params_0 = vcat( ones(sum(coordinates_free_params_B)), # loadings ones(1+n_trimmed), # variances of the cycles' innovations - 1e-4*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) + 1e-5*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ) # Lower bound @@ -335,7 +335,7 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e params_ub = vcat( +10*ones(sum(coordinates_free_params_B)), # loadings 1e+4*ones(1+n_trimmed), # variances of the cycles' innovations - 1e-2*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) + 1e-4*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); # Maximum likelihood From c6c89bc0e07d71df2957b662a3026159c90e5a62 Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Sat, 20 Jan 2024 17:16:32 +0000 Subject: [PATCH 13/14] Bump --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ae901f7..d733953 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MessyTimeSeriesOptim" uuid = "7115d47b-e118-4204-bc33-dbd7ed133e66" -version = "0.2.2" +version = "0.2.3" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" From f536faac7258944632f921e4049ad24dc4e07ffb Mon Sep 17 00:00:00 2001 From: fipelle <6272230+fipelle@users.noreply.github.com> Date: Tue, 13 Feb 2024 23:10:07 +0000 Subject: [PATCH 14/14] Numerically stable bounds --- src/initialisation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/initialisation.jl b/src/initialisation.jl index 0684643..d1d6c59 100644 --- a/src/initialisation.jl +++ b/src/initialisation.jl @@ -327,14 +327,14 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e # Lower bound params_lb = vcat( -10*ones(sum(coordinates_free_params_B)), # loadings - 1e-4*ones(1+n_trimmed), # variances of the cycles' innovations + 1e-2*ones(1+n_trimmed), # variances of the cycles' innovations 1e-6*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) ); # Upper bound params_ub = vcat( +10*ones(sum(coordinates_free_params_B)), # loadings - 1e+4*ones(1+n_trimmed), # variances of the cycles' innovations + 1e+2*ones(1+n_trimmed), # variances of the cycles' innovations 1e-4*ones(estim.n_trends) # variance of the trends' innovations (as a function of the cycles) );