Merge pull request #22 from fipelle/dev

Update to 0.2.1

Project.toml

@@ -1,6 +1,6 @@
 name = "MessyTimeSeriesOptim"
 uuid = "7115d47b-e118-4204-bc33-dbd7ed133e66"
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

src/initialisation.jl

@@ -285,7 +285,7 @@ function initial_univariate_decomposition(data::Union{FloatVector, JVector{Float
     cycle = smoothed_states_matrix[3, :];
 
     # Return output
-    return trend, drift_or_trend_lagged, cycle, status;
+    return trend, drift_or_trend_lagged, cycle, minimizer_bounded, status;
 end
 
 """
@@ -296,8 +296,8 @@ This function returns an initial estimate of the non-stationary and stationary c
 If `is_rw_trend` is false this function models the trend as in Kitagawa and Gersch (1996, ch. 8).
 """
 function initial_univariate_decomposition_kitagawa(data::Union{FloatVector, JVector{Float64}}, lags::Int64, ε::Float64, is_rw_trend::Bool; sigma_lb::Vector{Float64}=[1e-3; 1e3], sigma_ub::Vector{Float64}=[1e3; 1e4])
-    trend, _, cycle, _ = initial_univariate_decomposition(data, lags, ε, is_rw_trend, false, sigma_lb=sigma_lb, sigma_ub=sigma_ub);
-    return trend, cycle;
+    trend, _, cycle, minimizer_bounded, _ = initial_univariate_decomposition(data, lags, ε, is_rw_trend, false, sigma_lb=sigma_lb, sigma_ub=sigma_ub);
+    return trend, cycle, minimizer_bounded;
 end
 
 """
@@ -308,8 +308,8 @@ This function returns an initial estimate of the non-stationary and stationary c
 If `is_rw_trend` is false this function models the trend as a local linear trend.
 """
 function initial_univariate_decomposition_llt(data::Union{FloatVector, JVector{Float64}}, lags::Int64, ε::Float64, is_rw_trend::Bool; sigma_lb::Vector{Float64}=[1e-3; 1e3], sigma_ub::Vector{Float64}=[1e3; 1e4])
-    trend, drift, cycle, _ = initial_univariate_decomposition(data, lags, ε, is_rw_trend, true, sigma_lb=sigma_lb, sigma_ub=sigma_ub);
-    return trend, drift, cycle;
+    trend, drift, cycle, minimizer_bounded, _ = initial_univariate_decomposition(data, lags, ε, is_rw_trend, true, sigma_lb=sigma_lb, sigma_ub=sigma_ub);
+    return trend, drift, cycle, minimizer_bounded;
 end
 
 """
@@ -341,13 +341,16 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e
 
     # Compute individual trends
     trends = zeros(n_trimmed, T_trimmed);
-    cycles = zeros(n_trimmed, T_trimmed); 
+    cycles = zeros(n_trimmed, T_trimmed);
+    trends_variance = zeros(n_trimmed);
+
     for i=1:n_trimmed
         verb_message(estim.verb, "Initialisation > NLopt.LN_SBPLX, variable $(i)");
         drifts_selection_ids = findall(view(estim.trends_skeleton, i, :) .!= 0.0); # (i, :) is correct since it iterates series-wise
         if length(drifts_selection_ids) > 0
             drifts_selection_id = drifts_selection_ids[findmax(i -> estim.drifts_selection[i], drifts_selection_ids)[2]];
-            trends[i, :], cycles[i, :] = initial_univariate_decomposition_kitagawa(Y_trimmed[i, :], estim.lags, estim.ε, estim.drifts_selection[drifts_selection_id]==0);
+            trends[i, :], cycles[i, :], minimizer_bounded = initial_univariate_decomposition_kitagawa(Y_trimmed[i, :], estim.lags, estim.ε, estim.drifts_selection[drifts_selection_id]==0);
+            trends_variance[i] = minimizer_bounded[1];
         else
             cycles[i, :] .= Y_trimmed[i, :];
         end
@@ -358,7 +361,8 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e
 
     # Compute common trends. `common_trends` is equivalent to `trends` if there aren't common trends to compute.
     common_trends = zeros(estim.n_trends, T_trimmed);
-
+    common_trends_variance = zeros(estim.n_trends);
+
     # Looping over each trend in `common_trends`
     for i=1:estim.n_trends
 
@@ -369,8 +373,15 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e
         # Get selection of series onto which compute the current common trend
         coordinates_current_block = series_loading_on_current_trend[.!series_loading_on_further_trends];
 
-        # Unadjusted estimate (dividing by the coefficients in `estim.trends_skeleton` allows to appropriately weight each series)
+        # Unadjusted estimate of the trend (dividing by the coefficients in `estim.trends_skeleton` allows to appropriately weight each series)
         common_trends[i, :] = mean(trends[coordinates_current_block, :] ./ estim.trends_skeleton[coordinates_current_block, i], dims=1);
+
+        # Unadjusted estimate of the variances per series
+        unadjusted_trends_variance = trends_variance[coordinates_current_block];
+        unadjusted_trends_variance ./= sum(estim.trends_skeleton[coordinates_current_block, :] .!= 0.0, dims=2)[:];
+
+        # Adjusted estimate of the initial variance for the common trend (dividing by the coefficients in `estim.trends_skeleton` allows to appropriately weight each series)
+        common_trends_variance[i] = mean(unadjusted_trends_variance ./ estim.trends_skeleton[coordinates_current_block, i].^2);
 
         # Update `trends`
         for j in series_loading_on_current_trend
@@ -382,5 +393,5 @@ function initial_detrending(Y_untrimmed::Union{FloatMatrix, JMatrix{Float64}}, e
     detrended_data .-= estim.trends_skeleton*common_trends;
 
     # Return output
-    return common_trends, detrended_data;
+    return common_trends, common_trends_variance, detrended_data;
 end