From fa46b35c93575d5fafeaef5417b3bb8775735eee Mon Sep 17 00:00:00 2001 From: Michael Hatcher Date: Fri, 4 Oct 2024 14:39:23 +0100 Subject: [PATCH] Add files via upload --- PF_insert.m | 14 +++++--- P_matrix.m | 19 +++++----- Pos_Def.m | 33 +++++------------- Print.m | 12 +++++++ Select_solution.m | 2 +- Solutions_insert.m | 37 ++++++++++++++++++++ Solutions_insert_App_2.m | 47 +++++++++++++++++++++++++ Solutions_insert_App_2_loop.m | 65 +++++++++++++++++++++++++++++++++++ 8 files changed, 190 insertions(+), 39 deletions(-) create mode 100644 Print.m create mode 100644 Solutions_insert.m create mode 100644 Solutions_insert_App_2.m create mode 100644 Solutions_insert_App_2_loop.m diff --git a/PF_insert.m b/PF_insert.m index fa19bd8..3ef6262 100644 --- a/PF_insert.m +++ b/PF_insert.m @@ -1,6 +1,8 @@ %PF_insert. Guess and verify part of the algorithm. Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). %Any errors are my own. Last updated: 12/04/2023. +%Housekeeping +mstar = zeros(N_guess,1); X_stack = NaN(nvar,1); X_sol = NaN(nvar,T_sim-1,N_guess); X_sol_exc = NaN(nvar-nx,T_sim-1,N_guess); ind_sol = NaN(N_guess,T_sim-1); @@ -78,23 +80,27 @@ X_star = F*[X; X_e; X_lag] + G*e_vec(:,t) + H; X_star_store(t) = X_star; - if X_star >= X1_min && ind(t) == 1 || X_star <= X1_min && ind(t) == 0 + if X_star > X1_min && ind(t) == 1 || X_star <= X1_min && ind(t) == 0 %Can use X_star => X1_min && ind(t) == 1 for knife-edge case Verify(t) = 1; else break end - %if X_stack(1,t) == max(X1_min,X_star) does not work well due to numerical precision - %For accuracy, can use vpa(.) + %if X_stack(1,t) == max(X1_min,X_star) does not seem to work well end - if sum(Verify) == T_sim-1 && min(X_stack(1,T_guess:T_sim-1)) > X1_min + if sum(Verify) == T_sim-1 && min(X_stack(1,T_guess+1:T_sim-1)) > X1_min X_sol(:,:,m) = X_stack(:,:); X_sol_exc(:,:,m) = X_stack(1:nvar-nx,:); mstar(m) = m; ind_sol(m,:) = ind(1:T_sim-1)'; X_star_sol(m,:) = X_star_store(1:T_sim-1); + + if ~isempty(not_P) && not_P==0 % + break %Stop search if M is a P-matrix and a solution is found + end + end end diff --git a/P_matrix.m b/P_matrix.m index 1c8b007..bfa8a45 100644 --- a/P_matrix.m +++ b/P_matrix.m @@ -1,18 +1,19 @@ -%Script to determine if M is a P-matrix. %Uses ptest.m written by Michael -%Tsatsomeros (https://www.math.wsu.edu/faculty/tsat/matlab.html). The code can also be used with ptest3.m, by the same author. -%Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). Any errors are my own. Last updated: 12/04/2023. +%Script to determine if M is a P-matrix. %Uses ptest.m and ptest3.m written by Michael Tsatsomeros +%(see the MATLAB depot at https://www.math.wsu.edu/faculty/tsat/matlab.html). +%This script is written by Michael Hatcher (m.c.hatcher@soton.ac.uk). Any errors are my own. Last updated: 12/04/2023. run Pos_Def.m -if not_P == 0 && is_pd == 0 +if not_P ~=1 && is_pd == 0 + + %ptest_val = ptest(M); + ptest_val = ptest3(M); %faster version - if ptest(M)==0 + if ptest_val==0 not_P = 1; + else + not_P = 0; end - - %if ptest3(M)==0 - % not_P = 1; - %end end diff --git a/Pos_Def.m b/Pos_Def.m index 5ff7025..2d18309 100644 --- a/Pos_Def.m +++ b/Pos_Def.m @@ -1,35 +1,18 @@ %Check if M matrix is general positive definite --> P-matrix -%Check if simple necessary conditions for a P-matrix are violated %Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). -%Any errors are my own. Last updated: 12/04/2023. +%Any errors are my own. Last updated: 03/10/2024. -%Based on Holden (2022, Lemma 1 Part 2) +%Based on Holden (2023, Lemma 1 Part 2) eig_vals = eig(M + M'); length_M = length(M); -is_pd = nnz(eig_vals>0) == length_M; % flag to check if it is PD +is_pd = nnz(eig_vals>0) == length_M; % flag to check if it is PD, nnz - no. of non-zero elements -if is_pd == 0 +if is_pd == 1 - if det(M) <= 0 || min(diag(M)) <=0 - not_P = 1; + not_P = 0; - else +elseif min(diag(M)) <=0 || det(M) <= 0 -n_skip = max(round(length_M/10),1); -max_iter = max(length(M)-1,1); - -for iter = 1:n_skip:max_iter - -det_M = det(M(1:iter,1:iter)); - - if det_M <= 0 - not_P = 1; - break - end - -end - - end - -end + not_P = 1; +end \ No newline at end of file diff --git a/Print.m b/Print.m new file mode 100644 index 0000000..f5bff66 --- /dev/null +++ b/Print.m @@ -0,0 +1,12 @@ +%Print.m +%Messages to print to screen after using the algorithm + +if not_P == 0; disp('M matrix is a P-matrix'); +elseif not_P == 1; disp('M matrix is not a P-matrix'); +end + +if Message_1 == 1 + x_disp = sprintf('Increase N_guess if you want to include all %d guesses generated by Guesses_master or Guesses_master_2',d0); disp(x_disp) +elseif Message_2 == 1 + x_disp = sprintf('More guesses made than are generated by the Guesses_master file used: %d versus %d. The excess guesses are random guesses of 0s and 1s.',N_guess,d0); disp(x_disp) +end \ No newline at end of file diff --git a/Select_solution.m b/Select_solution.m index 09bc1e2..8823952 100644 --- a/Select_solution.m +++ b/Select_solution.m @@ -1,4 +1,4 @@ -%Select a solution (resolving the indeterminacy). Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). +%Select a solution (Remark 1, resolving the indeterminacy). Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). %Any errors are my own. Last updated: 12/04/2023. %Pick a solution (resolving indeterminacy) diff --git a/Solutions_insert.m b/Solutions_insert.m new file mode 100644 index 0000000..3df142c --- /dev/null +++ b/Solutions_insert.m @@ -0,0 +1,37 @@ +%Solutions_insert + +X_exog = []; +mstar(mstar==0) = []; + +if isempty(mstar) + + disp('No solution found. Check T_guess and N_guess are not too small.') + +elseif ~isempty(mstar) + + X_sol = X_sol(:,:,mstar); + X_exog = X_sol(nvar-nx+1:nvar,:,1); + X_sol_exc = X_sol_exc(:,:,mstar); + + solutions = reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2)); + ind_solutions = ind_sol(mstar,:); + + %Solution paths + x_fin = unique(solutions,'rows','stable'); + sol_ind = unique(ind_solutions,'rows','stable'); + + %No. of solutions + k = length(x_fin(:,1))/(nvar-nx); + + if k==1 + disp('Unique solution found') + else + disp('Multiple solutions found') + end + + %Final solutions + ind_fin = sol_ind; + x_fin = [x_fin; X_exog]; %Put exogenous last + + +end diff --git a/Solutions_insert_App_2.m b/Solutions_insert_App_2.m new file mode 100644 index 0000000..f8ddc13 --- /dev/null +++ b/Solutions_insert_App_2.m @@ -0,0 +1,47 @@ +%Solutions_insert_App_2 + +mstar(mstar==0) = []; + +if isempty(mstar) + + disp('No solution found. Check T_guess and N_guess are not too small.') + +elseif ~isempty(mstar) + + X_sol = X_sol(:,:,mstar); + X_exog = X_sol(nvar-nx+1:nvar,:,1); + X_sol_exc = X_sol_exc(:,:,mstar); + + %Note: The below is an ad-hoc workaround. Because in the model i(t) = i*(t) in the reference regime (slack), + %some of these paths will removed by the 'unique' function. To prevent this we add small positive number to + %final simulated point of variable 1 (i.e. i(t)). + for ii=1:length(X_sol(1,1,:)) + rng(1), noise = abs(randn)*1e-18; + orig = X_sol(1,T_sim-1,ii); orig_exc = X_sol_exc(1,T_sim-1,ii); + X_sol(1,T_sim-1,ii) = orig + noise/4; + X_sol_exc(1,T_sim-1,ii) = orig_exc + noise/4; + end + + solutions = reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2)); + ind_solutions = ind_sol(mstar,:); + + %Solution paths + x_fin = unique(solutions,'rows','stable'); + sol_ind = unique(ind_solutions,'rows','stable'); + + %No. of solutions + k = length(x_fin(:,1))/(nvar-nx); + + if k==1 + disp('Unique solution found') + else + disp('Multiple solutions found') + end + + %Final solutions + ind_fin = sol_ind; + x_fin = [x_fin; X_exog]; %Put exogenous last + +end + + diff --git a/Solutions_insert_App_2_loop.m b/Solutions_insert_App_2_loop.m new file mode 100644 index 0000000..dcd9d70 --- /dev/null +++ b/Solutions_insert_App_2_loop.m @@ -0,0 +1,65 @@ +%Solutions_insert_App_2_loop + +mstar(mstar==0) = []; + +if isempty(mstar) + + disp('No solution found. Check T_guess and N_guess are not too small.') + count_no_sol(j) = 1; + count_unique(j) = 0; + +elseif ~isempty(mstar) + + X_sol = X_sol(:,:,mstar); + X_exog = X_sol(nvar-nx+1:nvar,:,1); %Order exog last in x vector + X_sol_exc = X_sol_exc(:,:,mstar); + X_star_sol = X_star_sol(mstar,:); + + %Note: The below is an ad-hoc workaround. Because in the model i(t) = i*(t) in the reference regime (slack), + %some of these paths will removed by the 'unique' function. To prevent this we add small positive number to + %final simulated point of variable 1 (i.e. i(t)). + for ii=1:length(X_sol(1,1,:)) + rng(1), noise = abs(randn)*1e-18; + orig = X_sol(1,T_sim-1,ii); orig_exc = X_sol_exc(1,T_sim-1,ii); + X_sol(1,T_sim-1,ii) = orig + noise/4; + X_sol_exc(1,T_sim-1,ii) = orig_exc + noise/4; + end + + solutions = reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2)); + ind_solutions = ind_sol(mstar,:); + + %Solution paths + x_fin = unique(solutions,'rows','stable'); + sol_ind = unique(ind_solutions,'rows','stable'); + + count_no_sol(j) = 0; + + %No. of solutions + k = length(x_fin(:,1))/(nvar-nx); + + %Solutions for sim j + X_star_fin = unique(X_star_sol,'rows','stable'); + ind_fin = sol_ind; + x_fin = [x_fin; X_exog]; %Put exogenous last + + if k==1 + %disp('Unique solution found') + count_multi(j) = 0; count_unique(j) = 1; + count(j) = sum(ind_fin(1,:)==0); + + elseif k>=2 + %disp('Warning: Multiple solutions or no solution') + count_multi(j) = 1; + if k==2 + count(j) = sum(ind_fin(2,:)==0); %Count periods at bound + else + disp('Warning: More than two solutions') + end + + end + +end + + + +