new ; cls; closeall ; // -------------------------------------------------------- // // mequidata_endpaper_example_joint_identification_logit_14112018.gss // // This program considers the 2 x 2 x 2 model for the // example at the end of the paper. // The purpose of this example is showing that there is // a continuum of DGPs with joint identification but // without sequential identification. // // 2x2x2 means two-plaer, binary choice, two brances of unobserved heterogeneity // // The distribution of the private information is Logistic // // "Identification of Games of Incomplete Information with // Multiple Equilibria and Common Unobserved Heterogeneity," // by Aguirregabiria and Mira (QE) // // by Victor Aguirregabiria // // First version: November 2, 2018 // This version: November 14, 2018 // // -------------------------------------------------------- library pgraph gauss; wdir = "c:\\Users\\Victor\\Dropbox\\MYPAPERS\\GENETIC\\progau\\endpaper_example_joint_identification\\" ; buff = changedir(wdir) ; fileout = "mequidata_endpaper_example_joint_identification_logit_14112018.out" ; output file = ^fileout reset ; format /mb1 /ros 16,4 ; // --------------------------------------------------------------------- // MODEL // // - Define the CDF of the Logistic as: // cdflogit(x) = exp(x)/(1+exp(x)) // And the quantile function for the logistic distribution as: // qlogit(p) = ln(p) - ln(1-p) // // - For a given value of (z1,z2,kappa) the equilibrium equations are: // p1 = cdflogit(alpha_1 + beta_1 * p2) // p2 = cdflogit(alpha_2 + beta_2 * p1) // where p's are the CCPs, and alphas and betas are the payoffs. // // - We can write this system also as: // qlogit(p1) = alpha_1 + beta_1 * p2 // qlogit(p2) = alpha_2 + beta_2 * p1 // // - Therefore, for given alphas and betas, we can solve for the equilibrium CCPs // // - Each zi variable takes numz values, and kappa two values. // Therefore, there are numz * numz * 2 types of markets. // // --------------------------------------------------------------------- // -------------------------------- // 1. Parameters and constants // -------------------------------- nplayer = 2; // number of players numz = 4; // Number values for zi numkappa = 2; // Number of values of kappa zval = seqa(0,1,numz)/(numz-1); // Support of zi kappaval = seqa(0,1,numkappa)/(numkappa-1); // Support of kappa nmarket = numz*numz*numkappa; // Number of market types // {hA_mat, newseed} = rndu(numz,numz,5333799); // matrix with hA(z1,z2) = Pr(kappa = A | z1,z2) hA_mat = 0.7*ones(numz,numz); // matrix with hA(z1,z2) = Pr(kappa = A | z1,z2) mat_a1a2 = ( 0 ~ 0 ) | ( 0 ~ 1 ) | ( 1 ~ 0 ) | ( 1 ~ 1 ); numeqiter = 1000; // Maximum number of iterations to compute equilibrium // Payoff functions sigmapi = 1; a1_0 = -2.00; a1_z = 5.00; a1_k = 1.00; b1_0 = -2.00; b1_z = 0.00; b1_k = 0.00; a2_0 = -3.00; a2_z = 4.00; a2_k = 0.50; b2_0 = -3.00; b2_z = 0.00; b2_k = 0.00; alpha_1 = a1_0 + a1_z * zval + a1_k * (kappaval'); beta_1 = b1_0 + b1_z * zval + b1_k * (kappaval'); alpha_2 = a2_0 + a2_z * zval + a2_k * (kappaval'); beta_2 = b2_0 + b2_z * zval + b2_k * (kappaval'); /* // Random draws of payoff values {alpha_1, newseed} = rndn(numz,numkappa,3100); alpha_1 = sigmapi * alpha_1; {beta_1, newseed} = rndn(numz,numkappa,3200); beta_1 = sigmapi * beta_1; {alpha_2, newseed} = rndn(numz,numkappa,3300); alpha_2 = sigmapi * alpha_2; {beta_2, newseed} = rndn(numz,numkappa,3400); beta_2 = sigmapi * beta_2; */ // --------------------------------------------- // 2. Computing equilibrium probabilities // --------------------------------------------- // Matrix of equilibrium probabilities: // pequil_A for unobserved type A // pequil_B for unobserved type B // Column 1: equil prob for player 1 // Column 2: equil prob for player 2 // The rows are the different observed market types, // sorted by z1, and z2, in this order pequil_A = zeros(numz*numz,nplayer); pequil_B = zeros(numz*numz,nplayer); z1=1; do while z1<=numz; z2=1; do while z2<=numz; conta = (z1-1)*numz + z2; // Type kappa = A k=1; a1 = alpha_1[z1,k]; b1 = beta_1[z1,k]; a2 = alpha_2[z2,k]; b2 = beta_2[z2,k]; p_initial = zeros(2,1); p_final = zeros(2,1); eqiter=1; do while eqiter<=numeqiter; p_final = (a1 + b1 * p_initial[2]) | (a2 + b2 * p_initial[1]); p_final = exp(p_final)./(1+exp(p_final)); p_initial = p_final; eqiter = eqiter+1; endo; pequil_A[conta,.] = p_final'; // Type kappa = B k=2; a1 = alpha_1[z1,k]; b1 = beta_1[z1,k]; a2 = alpha_2[z2,k]; b2 = beta_2[z2,k]; p_initial = zeros(2,1); p_final = zeros(2,1); eqiter=1; do while eqiter<=numeqiter; p_final = (a1 + b1 * p_initial[2]) | (a2 + b2 * p_initial[1]); p_final = exp(p_final)./(1+exp(p_final)); p_initial = p_final; eqiter = eqiter+1; endo; pequil_B[conta,.] = p_final'; z2=z2+1; endo; z1=z1+1; endo; // Figure of equilibrium probabilities graphset ; _plwidth = 15 ; _pdate = "" ; _pcolor = {0,0}; xlabel("INDEX FOR STATE (z1,z2)") ; ylabel("EQUILIBRIUM PROBABILITY PLAYER 1") ; _plegctl = {1 5 10 0.1 } ; _plegstr = "Kappa = 0 \000Kappa = 1" ; xy(seqa(1,1,numz*numz),pequil_A[.,1]~pequil_B[.,1]); ylabel("EQUILIBRIUM PROBABILITY PLAYER 2") ; _plegctl = {1 5 10 0.5 } ; xy(seqa(1,1,numz*numz),pequil_A[.,2]~pequil_B[.,2]); // ----------------------------------------------------------- // 3. Jacobian matrix for joint identification // ----------------------------------------------------------- // The Jacobian matrix for joint identification has the following structure: // J_joint = [ Hlike ~ 0 ] // [ Dc_hP ~ Dc_pi ] // where: // - Hlike is the Hessian of the log-likelihod, with respect to h and P // - Dc_hP is the Jacobian of the constraints, with respect to h and P // - Dc_pi is the Jacobian of the constraints with respect to pi // // The parameters are sorted by (z1,z2), and then for each value of // (z1,z2) the are sorted as (hA,p1A,p1B,p2A,p2B) // ----------------------------------------------------------- // ----------------------------------------------------------- // 3.1. Obtaining the true distribution of Q(a1,a2|z1,z2) // ----------------------------------------------------------- // rows: values of (z1,z2); colsL values of (a1,a2) // The cols are sorted as: (0,0), (0,1), (1,0), (1,1) trueQ = zeros(numz*numz,4); z1=1; do while z1<=numz; z2=1; do while z2<=numz; conta = (z1-1)*numz + z2; p1_A = pequil_A[conta,1]; p2_A = pequil_A[conta,2]; p1_B = pequil_B[conta,1]; p2_B = pequil_B[conta,2]; hA = hA_mat[z1,z2]; trueQ[conta,1] = hA*(1-p1_A)*(1-p2_A) + (1-hA)*(1-p1_B)*(1-p2_B); trueQ[conta,2] = hA*(1-p1_A)*p2_A + (1-hA)*(1-p1_B)*p2_B; trueQ[conta,3] = hA*p1_A*(1-p2_A) + (1-hA)*p1_B*(1-p2_B); trueQ[conta,4] = hA*p1_A*p2_A + (1-hA)*p1_B*p2_B; z2=z2+1; endo; z1=z1+1; endo; xy(seqa(1,1,numz*numz),trueQ); // ---------------------------------------------------------------------- // 3.2. Obtaining the true scores dQ(a1,a2|z1,z2)/d(hA,p1A,p1B,p2A,p2B) // ---------------------------------------------------------------------- // The parameters, and therefore the columns of the score are sorted // as (hA,p1A,p1B,p2A,p2B) score_00 = zeros(numz*numz,5); score_01 = zeros(numz*numz,5); score_10 = zeros(numz*numz,5); score_11 = zeros(numz*numz,5); z1=1; do while z1<=numz; z2=1; do while z2<=numz; conta = (z1-1)*numz + z2; p1_A = pequil_A[conta,1]; p2_A = pequil_A[conta,2]; p1_B = pequil_B[conta,1]; p2_B = pequil_B[conta,2]; hA = hA_mat[z1,z2]; // (a1,a2) = (0,0) a1=0; a2=0; d1_A = (1-a1)*(1-p1_A) + a1*p1_A; d1_B = (1-a1)*(1-p1_B) + a1*p1_B; d2_A = (1-a2)*(1-p2_A) + a2*p2_A; d2_B = (1-a2)*(1-p2_B) + a2*p2_B; score_00[conta,.] = (d1_A * d2_A - d1_B * d2_B) ~ (hA*(2*a1-1)*d2_A) ~ ((1-hA)*(2*a1-1)*d2_B) ~ (hA*d1_A*(2*a2-1)) ~ ((1-hA)*d1_B*(2*a2-1)) ; // (a1,a2) = (0,1) a1=0; a2=1; d1_A = (1-a1)*(1-p1_A) + a1*p1_A; d1_B = (1-a1)*(1-p1_B) + a1*p1_B; d2_A = (1-a2)*(1-p2_A) + a2*p2_A; d2_B = (1-a2)*(1-p2_B) + a2*p2_B; score_01[conta,.] = (d1_A * d2_A - d1_B * d2_B) ~ (hA*(2*a1-1)*d2_A) ~ ((1-hA)*(2*a1-1)*d2_B) ~ (hA*d1_A*(2*a2-1)) ~ ((1-hA)*d1_B*(2*a2-1)) ; // (a1,a2) = (1,0) a1=1; a2=0; d1_A = (1-a1)*(1-p1_A) + a1*p1_A; d1_B = (1-a1)*(1-p1_B) + a1*p1_B; d2_A = (1-a2)*(1-p2_A) + a2*p2_A; d2_B = (1-a2)*(1-p2_B) + a2*p2_B; score_10[conta,.] = (d1_A * d2_A - d1_B * d2_B) ~ (hA*(2*a1-1)*d2_A) ~ ((1-hA)*(2*a1-1)*d2_B) ~ (hA*d1_A*(2*a2-1)) ~ ((1-hA)*d1_B*(2*a2-1)) ; // (a1,a2) = (1,1) a1=1; a2=1; d1_A = (1-a1)*(1-p1_A) + a1*p1_A; d1_B = (1-a1)*(1-p1_B) + a1*p1_B; d2_A = (1-a2)*(1-p2_A) + a2*p2_A; d2_B = (1-a2)*(1-p2_B) + a2*p2_B; score_11[conta,.] = (d1_A * d2_A - d1_B * d2_B) ~ (hA*(2*a1-1)*d2_A) ~ ((1-hA)*(2*a1-1)*d2_B) ~ (hA*d1_A*(2*a2-1)) ~ ((1-hA)*d1_B*(2*a2-1)) ; z2=z2+1; endo; z1=z1+1; endo; // ---------------------------------------------------------------------- // 3.3. Checking that the sum of scores is zero // ---------------------------------------------------------------------- sumscore = score_00 + score_01 + score_10 + score_11; sumscore = sumc(sumscore)/rows(sumscore); ""; "Checking the gradient of the log-likelihood is zero"; "Gradient =";; sumscore'; ""; // ---------------------------------------------------------------------- // 3.4. Hessian of log-likelihood // ---------------------------------------------------------------------- // Hessian matrix: Block diagonal, where each block corresponds to a value of (z1,z2) // ∇²ℓ = ∑_{(a₁,a₂)∈{0,1}²} (1/(Q(a₁,a₂))) (∂Q(a₁,a₂)/∂[hA,P]) (∂Q(a₁,a₂)/∂[hA,P]′) // ---------------------------------------------------------------------- hessian_like = zeros(5*numz*numz,5*numz*numz); z1=1; do while z1<=numz; z2=1; do while z2<=numz; conta = (z1-1)*numz + z2; ind1 = (conta-1)*5 + 1; ind2 = conta*5; hessian_like[ind1:ind2,ind1:ind2] = (score_00[conta,.]'*score_00[conta,.])./trueQ[conta,1] + (score_01[conta,.]'*score_01[conta,.])./trueQ[conta,2] + (score_10[conta,.]'*score_10[conta,.])./trueQ[conta,3] + (score_11[conta,.]'*score_11[conta,.])./trueQ[conta,4]; z2=z2+1; endo; z1=z1+1; endo; // ---------------------------------------------------------------------- // 3.5. Eigenvalues of Hessian matrix // ---------------------------------------------------------------------- hess_eig = eigh(hessian_like); xy(seqa(1,1,rows(hess_eig)),hess_eig); rank_hess = rank(hessian_like); "" "Rank of Hessian of likelihood:";; rank_hess; ""; // ---------------------------------------------------------------------- // 3.6. Restrictions: c(i,z1,z2,kappa; pi, P) // Jacobian with respect to (h,P): Dc_hp // ---------------------------------------------------------------------- // - The parameters h and P are sorted in the same way as in the Hessian, // such that we can concatenate vertically the Hessian matrix and this Jacobian matrix. // - That is, the parameters are sorted by (z1,z2), and then for each value of // (z1,z2) they are sorted as (hA,p1A,p1B,p2A,p2B) // - There are as many constraints as parameters in P. // - The constraints are sorted first by the values of (z1,z2), // and then for each value of (z1,z2) they are sorted as (c1A,c1B,c2A,c2B) // - For each block, the system of four constraints is: // c1A: alpha_1A + beta_1A * p2A - qlogit(p1A) // c1B: alpha_1B + beta_1B * p2B - qlogit(p1B) // c2A: alpha_2A + beta_2A * p1A - qlogit(p2A) // c2B: alpha_2B + beta_2B * p1B - qlogit(p2B) // - Note that the column associated to parameter hA is zero because // these constraints do not depend on hA. // - Remember that qlogit(p) = ln(p) - ln(1-p). Therefore, the dertivative // of this function is 1/(p*(1-p)) // - Matrix Dc_hp is block diagonal. For each value of (z1,z2) we have // a "block" matrix of dimension 4 x 5 // - And the jacobian matrix is: // ( 0 ~ -1/(p1A*(1-p1A)) ~ 0 ~ beta_1A ~ 0 ) // | ( 0 ~ 0 ~ -1/(p1B*(1-p1B)) ~ 0 ~ beta_1B ) // | ( 0 ~ beta_2A ~ 0 ~ -1/(p2A*(1-p2A)) ~ 0 ) // | ( 0 ~ 0 ~ beta_2B ~ 0 ~ -1/(p2B*(1-p2B) ) // ------------------------------------------------------------------------------------------------------- Dc_hp = zeros(numz*numz*2*2,numz*numz*5); z1=1; do while z1<=numz; z2=1; do while z2<=numz; conta = (z1-1)*numz + z2; row_ind1 = (conta-1)*4 + 1; row_ind2 = conta*4; col_ind1 = (conta-1)*5 + 1; col_ind2 = conta*5; p1_A = pequil_A[conta,1]; p2_A = pequil_A[conta,2]; p1_B = pequil_B[conta,1]; p2_B = pequil_B[conta,2]; b1_A = beta_1[z1,1]; b1_B = beta_1[z1,2]; b2_A = beta_2[z2,1]; b2_B = beta_2[z2,2]; blockmat = ( 0 ~ -1/(p1_A*(1-p1_A)) ~ 0 ~ b1_A ~ 0 ) | ( 0 ~ 0 ~ -1/(p1_B*(1-p1_B)) ~ 0 ~ b1_B ) | ( 0 ~ b2_A ~ 0 ~ -1/(p2_A*(1-p2_A)) ~ 0 ) | ( 0 ~ 0 ~ b2_B ~ 0 ~ -1/(p2_B*(1-p2_B)) ); Dc_hp[row_ind1:row_ind2,col_ind1:col_ind2] = blockmat; z2=z2+1; endo; z1=z1+1; endo; // ---------------------------------------------------------------------- // 3.7. Restrictions: c(i,z1,z2,kappa; pi, P) // Rank of stacked matrix "Hessian | Dc_hp" // ---------------------------------------------------------------------- Hess_Dchp = hessian_like | Dc_hp; rank_hess_Dchp = rank(Hess_Dchp); ""; "Number of columns of Hessian =";; cols(hessian_like); " Rank of Hessian =";; rank(hessian_like); "Number of columns of Hess_Dchp =";; cols(Hess_Dchp); " Rank of Hess_Dchp =";; rank(Hess_Dchp); ""; // ---------------------------------------------------------------------- // 3.8. Restrictions: c(i,z1,z2,kappa; pi, P) // Jacobian with respect to pi: Dc_pi // ---------------------------------------------------------------------- // - The constraints c are sorted in the same way as in Jacobian matrix Dc_hp, // such that we can concatenate horizontally the matrices Dc_hp and Dc_pi. // - There are as many constraints as parameters in P. // - The constraints are sorted first by the values of (z1,z2), // and then for each value of (z1,z2) they are sorted as (c1A,c1B,c2A,c2B) // - For each block, the system of four constraints is: // c1A: alpha_1A + beta_1A * p2A - qlogit(p1A) // c1B: alpha_1B + beta_1B * p2B - qlogit(p1B) // c2A: alpha_2A + beta_2A * p1A - qlogit(p2A) // c2B: alpha_2B + beta_2B * p1B - qlogit(p2B) // - The payoff parameters, and so the columns of the Hacobian matrix Dc_pi, // are sorted first by player: all the payoff params of player 1, and then // all the payoff params of player 2. Then, for each player j, the are sorted // by the variable zj. Finally, for player j and a value zj, the columns // correspond to: (alpha_jA, beta_jA, alpha_jB, beta_jB) // // - The submatrix of Dc_pi associated player 1 and value z1, for rows (c1A, c1B, c2A, c2B) // and columns (alpha_1A, beta_1A, alpha_1B, beta_1B) is: // ( 1 ~ pequil_A[conta,2] ~ 0 ~ 0 ) // | ( 0 ~ 0 ~ 1 ~ pequil_B[conta,2] ) // | ( 0 ~ 0 ~ 0 ~ 0 ) // | ( 0 ~ 0 ~ 0 ~ 0 ) // - Similarly, the submatrix of Dc_pi associated player 2 and value z2, for rows (c1A, c1B, c2A, c2B) // and columns (alpha_2A, beta_2A, alpha_2B, beta_2B) is: // ( 0 ~ 0 ~ 0 ~ 0 ) // | ( 0 ~ 0 ~ 0 ~ 0 ) // | ( 1 ~ pequil_A[conta,1] ~ 0 ~ 0 ) // | ( 0 ~ 0 ~ 1 ~ pequil_B[conta,1] ) // ------------------------------------------------------------------------------------------------------- Dc_pi = zeros(numz*numz*2*2,numz*2*4); z1=1; do while z1<=numz; z2=1; do while z2<=numz; conta = (z1-1)*numz + z2; row_ind1 = (conta-1)*4 + 1; row_ind2 = conta*4; col_min_1 = (z1-1)*4 + 1; col_max_1 = z1*4; col_min_2 = numz*4 + (z2-1)*4 + 1; col_max_2 = numz*4 + z2*4; Dc_pi[row_ind1:row_ind2,col_min_1:col_max_1] = ( 1 ~ pequil_A[conta,2] ~ 0 ~ 0 ) | ( 0 ~ 0 ~ 1 ~ pequil_B[conta,2] ) | ( 0 ~ 0 ~ 0 ~ 0 ) | ( 0 ~ 0 ~ 0 ~ 0 ); Dc_pi[row_ind1:row_ind2,col_min_2:col_max_2] = ( 0 ~ 0 ~ 0 ~ 0 ) | ( 0 ~ 0 ~ 0 ~ 0 ) | ( 1 ~ pequil_A[conta,1] ~ 0 ~ 0 ) | ( 0 ~ 0 ~ 1 ~ pequil_B[conta,1] ); z2=z2+1; endo; z1=z1+1; endo; // ---------------------------------------------------------------------- // 4. Ranks of J_joint and Jseq and other Jacobian matrices // ---------------------------------------------------------------------- mat_zeros_right = zeros(rows(hessian_like),cols(Dc_pi)); mat_zeros_left = zeros(rows(Dc_pi),cols(hessian_like)); J_joint = (hessian_like ~ mat_zeros_right) | (Dc_hp ~ Dc_pi ); J_seq = (hessian_like ~ mat_zeros_right) | (mat_zeros_left ~ Dc_pi ); sel_hA = seqa(1,5,numz*numz); hessian_hA = hessian_like[.,sel_hA]; rank_hessian = rank(hessian_like); rank_hessian_hA = rank(hessian_hA); rank_Dchp = rank(Dc_hp); rank_hessian_Dchp = rank(hessian_like | Dc_hp); rank_Dcpi = rank(Dc_pi); rank_matzeros_Dcpi = rank(mat_zeros_right | Dc_pi); rank_jjoint = rank(J_joint); rank_jseq = rank(J_seq); ""; "----------------------------------------------------------------------"; "Rank of the Hessian"; "Number of columns of Hessian =";; cols(hessian_like); " Rank of Hessian =";; rank_hessian; "----------------------------------------------------------------------"; "Rank of the Hessian_hA"; "Number of columns of Hessian_hA =";; cols(hessian_hA); " Rank of Hessian_hA =";; rank_hessian_hA; "----------------------------------------------------------------------"; "Rank of the Jacobian Dc_hP"; "Number of columns of Jacobian Dc_hP =";; cols(Dc_hp); " Rank of Jacobian Dc_hP =";; rank_Dchp; "----------------------------------------------------------------------"; "Rank of the vertical concatenation fo Hessian and Jacobian Dc_hP"; "Number of columns of Hessian | Dc_hP =";; cols(Dc_hp); " Rank of Hessian | Dc_hP =";; rank_hessian_Dchp; "----------------------------------------------------------------------"; "Rank of the Jacobian Dc_pi"; "Number of columns of Jacobian Dc_pi =";; cols(Dc_pi); " Rank of Jacobian Dc_pi =";; rank_Dcpi; "----------------------------------------------------------------------"; "Rank of the vertical concatenation fo Mat zeros and Jacobian Dc_pi"; "Number of columns of Mat zeros | Dc_pi =";; cols(Dc_pi); " Rank of Mat zeros | Dc_pi =";; rank_matzeros_Dcpi; "----------------------------------------------------------------------"; "Rank of complete Jacobian matrix J_joint"; "Number of columns of J_joint =";; cols(J_joint); " Rank of J_joint =";; rank_jjoint; "----------------------------------------------------------------------"; "Rank of complete Jacobian matrix J_seq"; "Number of columns of J_seq =";; cols(J_seq); " Rank of J_seq =";; rank_jseq; "----------------------------------------------------------------------"; ""; // ---------------------------------------------------------------------- // 5. Condition numbers of J_joint' * J_joint and other matrices // ---------------------------------------------------------------------- // Condition number of J_joint' * J_joint JJjoint_mat = J_joint' * J_joint; eig_JJ = eig(JJjoint_mat); abs_eig_JJ = abs(eig_JJ); abs_eig_JJ = sortc(abs_eig_JJ,1); minabs_eig_JJjoint = abs_eig_JJ[1]; maxabs_eig_JJjoint = abs_eig_JJ[rows(abs_eig_JJ)]; condnum_JJjoint = maxabs_eig_JJjoint/minabs_eig_JJjoint; // Condition number of J_seq' * J_seq JJseq_mat = J_seq' * J_seq; eig_seq = eig(JJseq_mat); abs_eig_seq = abs(eig_seq); abs_eig_seq = sortc(abs_eig_seq,1); minabs_eig_seq = abs_eig_seq[1]; maxabs_eig_seq = abs_eig_seq[rows(abs_eig_seq)]; condnum_seq = maxabs_eig_seq/minabs_eig_seq; // Condition number of Dc_pi'* Dc_pi sqDc_pi = Dc_pi'*Dc_pi ; eig_sqDcpi = eig(sqDc_pi); abs_eig_sqDcpi = abs(eig_sqDcpi); abs_eig_sqDcpi = sortc(abs_eig_sqDcpi,1); minabs_eig_sqDcpi = abs_eig_sqDcpi[1]; maxabs_eig_sqDcpi = abs_eig_sqDcpi[rows(abs_eig_sqDcpi)]; condnum_sqDcpi = maxabs_eig_sqDcpi/minabs_eig_sqDcpi; // Condition number of HD'*HD // where HD matrix is the vertical concatenation Hessian | Dc_hp HD = hessian_like | Dc_hp; sqHD = HD'*HD; eig_sqHD = eig(sqHD); abs_eig_sqHD = abs(eig_sqHD); abs_eig_sqHD = sortc(abs_eig_sqHD,1); minabs_eig_sqHD = abs_eig_sqHD[1]; maxabs_eig_sqHD = abs_eig_sqHD[rows(abs_eig_sqHD)]; condnum_sqHD = maxabs_eig_sqHD/minabs_eig_sqHD; // Condition number of Hessian_h' * Hessian_h // where Hessian_h is the matrix with the columns // of the Hessian associated to hA sel_hA = seqa(1,5,numz*numz); hessian_hA = hessian_like[.,sel_hA]; sqhess_hA = hessian_hA'*hessian_hA; eig_sqhess_hA = eig(sqhess_hA); abs_eig_sqhess_hA = abs(eig_sqhess_hA); abs_eig_sqhess_hA = sortc(abs_eig_sqhess_hA,1); minabs_eig_sqhess_hA = abs_eig_sqhess_hA[1]; maxabs_eig_sqhess_hA = abs_eig_sqhess_hA[rows(abs_eig_sqhess_hA)]; condnum_sqhess_hA = maxabs_eig_sqhess_hA/minabs_eig_sqhess_hA; // Condition number of Dc_P' * Dc_P // where Dc_P is the Jacobian matrix associated to the // derivatives with respect the CCPs seq_1234 = ones(numz*numz,1) .*. (1|2|3|4); sel_p = seqa(1,5,numz*numz) .*. (1|1|1|1); sel_p = seq_1234 + sel_p; Dc_p = Dc_hp[.,sel_p]; sqDcp = Dc_p'*Dc_p; eig_sqDcp = eig(sqDcp); abs_eig_sqDcp = abs(eig_sqDcp); abs_eig_sqDcp = sortc(abs_eig_sqDcp,1); minabs_eig_sqDcp = abs_eig_sqDcp[1]; maxabs_eig_sqDcp = abs_eig_sqDcp[rows(abs_eig_sqDcp)]; condnum_sqDcp = maxabs_eig_sqDcp/minabs_eig_sqDcp; ""; "----------------------------------------------------------------------------"; "Condition number of Hessian_h' * Hessian_h"; " Minimum absolute eigenvalue =";; minabs_eig_sqhess_hA; " Maximum absolute eigenvalue =";; maxabs_eig_sqhess_hA; "Condition Number of Hessian_h' * Hessian_h =";; condnum_sqhess_hA; "----------------------------------------------------------------------------"; ""; "----------------------------------------------------------------------------"; "Condition number of Dc_p' * Dc_p"; " Minimum absolute eigenvalue =";; minabs_eig_sqDcp; " Maximum absolute eigenvalue =";; maxabs_eig_sqDcp; "Condition Number of Dc_p' * Dc_p =";; condnum_sqDcp; "----------------------------------------------------------------------------"; ""; "----------------------------------------------------------------------------"; "Condition number of Dc_pi' * Dc_pi"; " Minimum absolute eigenvalue =";; minabs_eig_sqDcpi; " Maximum absolute eigenvalue =";; maxabs_eig_sqDcpi; " Condition Number of Dc_pi' * Dc_pi =";; condnum_sqDcpi; "----------------------------------------------------------------------------"; ""; "----------------------------------------------------------------------------"; "Condition number of HD' * HD, where HD = Hessian | Dc_hp"; " Minimum absolute eigenvalue =";; minabs_eig_sqHD; " Maximum absolute eigenvalue =";; maxabs_eig_sqHD; " Condition Number of HD' * HD =";; condnum_sqHD; "----------------------------------------------------------------------------"; ""; "----------------------------------------------------------------------------"; "Condition number of J_joint' * J_joint"; " Minimum absolute eigenvalue =";; minabs_eig_JJjoint; " Maximum absolute eigenvalue =";; maxabs_eig_JJjoint; "Condition Number of J_joint' * J_joint =";; condnum_JJjoint; "----------------------------------------------------------------------------"; ""; "----------------------------------------------------------------------------"; "Condition number of J_seq' * J_seq"; " Minimum absolute eigenvalue =";; minabs_eig_seq; " Maximum absolute eigenvalue =";; maxabs_eig_seq; "Condition Number of J_joint' * J_joint =";; condnum_seq; "----------------------------------------------------------------------------"; ""; numauto = 10; xy(seqa(1,1,numauto),abs_eig_JJ[1:numauto]); closeall; end ;