1. stevejb
  2. is-solver

Commits

stevejb  committed 0d8e3c3

Some optimizations so that DDW smm is absurdly fast

  • Participants
  • Parent commits ae32c88
  • Branches default

Comments (0)

Files changed (1)

File solve_isC.cpp

View file
  • Ignore whitespace
     controlTrajectory(n,0) = iRand;
   }
 
+#pragma omp parallel for  
   for(int n = 0; n < Nsim; n++) {
     for(int t = 1; t < T; t++) {
       controlTrajectory(n,t) = policy(shockTrajectory(n,t), controlTrajectory(n, t-1));
 
-      if(DBGPRINT) {      
-      // to shed some light
-      if( (n == 0 ) && (t > 0) && (t < 200) ) {
-        cout << "t: " << t <<  ", SHOCK:" << shockTrajectory(n,t) << ", CTRL PREV: " << controlTrajectory(n, t-1) << endl;
-        cout << "     CTRL CHOICE: " <<       controlTrajectory(n,t) << endl;
-      }
-      }
+      // if(DBGPRINT) {      
+      // // to shed some light
+      // if( (n == 0 ) && (t > 0) && (t < 200) ) {
+      //   cout << "t: " << t <<  ", SHOCK:" << shockTrajectory(n,t) << ", CTRL PREV: " << controlTrajectory(n, t-1) << endl;
+      //   cout << "     CTRL CHOICE: " <<       controlTrajectory(n,t) << endl;
+      // }
+      // }
 
     }
-  }
+  } // end for n
   return(controlTrajectory.middleCols(T-T_prev, T_prev));
 
 }
 //   printSimSummary(sim_control, pathsim, controlstates, shockstates,
 //                   fname, par, value, NULL);
 // }
-void printSimSummary(MatrixXi sim_control, MatrixXi pathsim,
+MapStringDoubleType printSimSummary(MatrixXi sim_control, MatrixXi pathsim,
                      MatrixRMXd controlstates, MatrixRMXd shockstates,
                      string fname, Params par,
                      MatrixRMXd value,
   mf2.close();
 
 
+  MapStringDoubleType tr_moments_map;
+  tr_moments_map["mean_inv"]                = mean_inv;
+  tr_moments_map["var_inv"]                 = var_inv;
+  tr_moments_map["mean_lev"]                = mean_lev;
+  tr_moments_map["var_lev"]                 = var_lev;
+  tr_moments_map["mean_equity_issuance"]    = mean_equity_issuance;
+  tr_moments_map["var_equity_issuance"]     = var_equity_issuance;
+  tr_moments_map["mean_tobins_q"]           = mean_tobins_q;
+  tr_moments_map["var_tobins_q"]            = var_tobins_q;
+  tr_moments_map["mean_cash_bal"]           = mean_cash_bal;
+  tr_moments_map["var_cash_bal"]            = var_cash_bal;
+  tr_moments_map["mean_opr_income"]         = mean_opr_income;
+  tr_moments_map["var_opr_income"]          = var_opr_income;
+  tr_moments_map["ser_cor"]                 = opr_income_reg.beta;
+  tr_moments_map["var_opr_income_innov"]    = opr_income_reg.var_resid;
+  tr_moments_map["inv_3rd_central"]         = inv_3rd_central;
+  tr_moments_map["inv_skewness"]            = inv_skewness;
+  tr_moments_map["inv_kurtosis"]            = inv_kurtosis;
+  tr_moments_map["lev_3rd_central"]         = lev_3rd_central;
+  tr_moments_map["lev_skewness"]            = lev_skewness;
+  tr_moments_map["lev_kurtosis"]            = lev_kurtosis;
+  tr_moments_map["innov_inc_3rd_central"]   = innov_inc_3rd_central;
+  tr_moments_map["innov_inc_skewness"]      = innov_inc_skewness;
+  tr_moments_map["innov_inc_kurtosis"]      = innov_inc_kurtosis;
+  tr_moments_map["opr_inc_3rd_central"]     = opr_inc_3rd_central;
+  tr_moments_map["opr_inc_skewness"]        = opr_inc_skewness;
+  tr_moments_map["opr_inc_kurtosis"]        = opr_inc_kurtosis;
+  tr_moments_map["cash_bal_3rd_central"]    = cash_bal_3rd_central;
+  tr_moments_map["cash_bal_skewness"]       = cash_bal_skewness;
+  tr_moments_map["cash_bal_kurtosis"]       = cash_bal_kurtosis;
+  tr_moments_map["tobinq_q_3rd_central"]    = tobinq_q_3rd_central;
+  tr_moments_map["tobinq_q_skewness"]       = tobinq_q_skewness;
+  tr_moments_map["tobinq_q_kurtosis"]       = tobinq_q_kurtosis;
+  tr_moments_map["eq_issuance_3rd_central"] = eq_issuance_3rd_central;
+  tr_moments_map["eq_issuance_skewness"]    = eq_issuance_skewness;
+  tr_moments_map["eq_issuance_kurtosis"]    = eq_issuance_kurtosis;
+  tr_moments_map["cov_inv_auto1"]           = cov_inv_auto1;
+  tr_moments_map["cov_inv_auto2"]           = cov_inv_auto2;
+  tr_moments_map["cov_inv_auto3"]           = cov_inv_auto3;
+  tr_moments_map["cov_opinc_auto1"]         = cov_opinc_auto1;
+  tr_moments_map["cov_opinc_auto2"]         = cov_opinc_auto2;
+  tr_moments_map["cov_opinc_auto3"]         = cov_opinc_auto3;
+  tr_moments_map["cov_inv_opinc"]           = cov_inv_opinc;
+  tr_moments_map["cov_kl3_opinc3"]          = cov_kl3_opinc3;
+  tr_moments_map["cov_lev_auto1"]           = cov_lev_auto1;
+  tr_moments_map["cov_lev_auto2"]           = cov_lev_auto2;
+  tr_moments_map["cov_lev_auto3"]           = cov_lev_auto3;
+  tr_moments_map["cov_lev_opincL1"]         = cov_lev_opincL1;
+  tr_moments_map["cov_lev_opincL2"]         = cov_lev_opincL2;
+  tr_moments_map["cov_lev_opincL3"]         = cov_lev_opincL3;
+  tr_moments_map["soa_beta_icept"]          = beta_soa(0,0);
+  tr_moments_map["soa_beta_lev1"]           = beta_soa(1,0);
+  tr_moments_map["soa_beta_val"]            = beta_soa(2,0);
+  tr_moments_map["soa_beta_cf"]             = beta_soa(3,0);
+
+  return(tr_moments_map);
   // printDecMatrix(operating_income_MOMENT, fname_base + "oper_income.csv");
   // printDecMatrix(opr_income_reg.residuals, fname_base + "oper_income_resids.csv");
   // printDecMatrix(eq_issuance_MOMENT, fname_base + "eq_issuance_MOMENT.csv");
 					       par, V_new, mapping, fname);
   } else {
     cout << "PSS No Trend " << endl;
-    printSimSummary(sim_control, pathsim, controlstates, shockstates, fname + "full-sim-notrend.txt", par, V_new, mapping, fname);
+    simulatedMoments = printSimSummary(sim_control, pathsim, controlstates, shockstates, fname + "full-sim-notrend.txt", par, V_new, mapping, fname);
   }
   cout << "...done. See: " << fname + "full-sim.txt" << endl << endl;
   cout << "========================================" << endl;
 
-  if(trended) {
-    return(simulatedMoments);
-  }
+  return(simulatedMoments);
+  // if(trended) {
+  //   return(simulatedMoments);
+  // }
   
-  MapStringDoubleType badmap;
-  badmap["NULL"] = 1;
-  return(badmap);
+  // MapStringDoubleType badmap;
+  // badmap["NULL"] = 1;
+  // return(badmap);
 
 } // end simulateAndMoments
 
   // fn_sub_mat is the initial random guess of the Bellman equation
   MatrixRMXd staticmat = MatrixRMXd(shockstates.rows()*statmatlen, statmatlen);
 
-  cout << "do bellman" << endl;
-  #pragma omp parallel for
-  for(int i = 0; i < shockstates.rows(); i++) {
-    cout << "iter: " << i << ", " << shockstates(i,0) << ", " << shockstates(i,1) << endl;
-    int outerrow = statmatlen * i;
-    staticmat.block(outerrow, 0, statmatlen, statmatlen) =
-      bellmanEq(par, k, k_prime, p, p_prime, shockstates(i,0), shockstates(i,1), i, SAVE_ALL, fname);
+  if(!(rtp.do_smm_for_ddw)) {
+    cout << "do bellman" << endl;
+#pragma omp parallel for
+    for(int i = 0; i < shockstates.rows(); i++) {
+      cout << "iter: " << i << ", " << shockstates(i,0) << ", " << shockstates(i,1) << endl;
+      int outerrow = statmatlen * i;
+      staticmat.block(outerrow, 0, statmatlen, statmatlen) =
+	bellmanEq(par, k, k_prime, p, p_prime, shockstates(i,0), shockstates(i,1), i, SAVE_ALL, fname);
+    }
+    cout << "done bellman" << endl;
+    cout << "staticmat: " << staticmat.rows() << " x " << staticmat.cols() << endl;
   }
-  cout << "done bellman" << endl;
-  cout << "staticmat: " << staticmat.rows() << " x " << staticmat.cols() << endl;
-
 
   // GENERATE THE RANDOM GUESS
   // cout << "start random" << endl;
   ///////////////////////////////////////////////////////////////
   // THIS IS FOR THE MAIN PROBLEM
   string fname2 = fname;
-  bool cvg_success;
-  // POPULATE THE INITIAL VALUE FUNCTION
-  valueThreaded(Nshocks, nstates, V_new, fullmat);
-  cvg_success = vfiter(T_BIG, staticmat, fullmat, V_new, policy, fname, par, iterpar);
-  printDecMatrix(V_new, fname + "value-barr.csv");
-
+  bool cvg_success = false;
+
+  if(!(rtp.do_smm_for_ddw)) {
+    // POPULATE THE INITIAL VALUE FUNCTION
+    valueThreaded(Nshocks, nstates, V_new, fullmat);
+    cvg_success = vfiter(T_BIG, staticmat, fullmat, V_new, policy, fname, par, iterpar);
+    printDecMatrix(V_new, fname + "value-barr.csv");
+  }
   if(!(rtp.do_smm)) {
     fname2 = fname + "value-barr.bin";
     WriteMatrix(fname2.c_str(), V_new);
   momentKeys = initializeMomentKeys();
 
 
-  if(cvg_success) {
+  if(cvg_success & !(rtp.do_smm_for_ddw)) {
 
     ////////////////////////////////////////
     // FULL POLICY SIM
   // cout << "========================================" << endl << endl;
   // cout << "Out early!" << endl;
 
-
-  if( (rtp.do_trans_only == true) && ( cvg_success == true) ) {
+  MapStringDoubleType ddw_moment_map;
+  if( ( (rtp.do_trans_only == true) && ( cvg_success == true) ) ||
+      rtp.do_smm_for_ddw == true ) {
     ////////////////////////////////////////////////////////////
     // NOW, SOLVE A SUBSET OF THE MAIN PROBLEM
     // FOR THE TRANSIENT ONLY CASE
 
 
     srand(101);
+    if(rtp.do_smm_for_ddw) {
+      cout << "Zoomed right to DDW estimation" << endl;
+    }
 
     // meaning the vector Z
 
       Nsim = 20000;
       T = 30;
       cout << "Nsim: " << Nsim << ", T: " << T << endl;
-      simulateAndMoments( Zprob, Nsim, T, fname_trans, policy_trans, initial_low_k_p,
+      ddw_moment_map = simulateAndMoments( Zprob, Nsim, T, fname_trans, policy_trans, initial_low_k_p,
                           par, controlstates, shockstates, V_trans, transonly_mapping);
 
     } else {
     // PARAMETER rtp.momentlist
     cout << "MOMENTS:" << endl;
     for(int iiii = 0; iiii < ddw_smm_Nmoments; iiii++) { 
-      ddw_smm_momentmat(iiii,0) = trendedMoments[ddw_smm_momentvec[iiii]];
+      ddw_smm_momentmat(iiii,0) = ddw_moment_map[ddw_smm_momentvec[iiii]];
       cout << setw(3) <<  iiii  << ":" << setw(10) <<  ddw_smm_momentvec[iiii] << " : " 
 	   << setw(8) << ddw_smm_momentmat(iiii,0) << endl;
     }