C++ code for gsl minsky simulation

Below follows the c++ code for the minsky simulation program.

// ~/xlin/gsl/testminsky.cc
// 150924 003 Task now -> grid
// Basert på 150921 001 testminskyall.cc
// Unified source for (74) and (77) variants
//
// Todo:  Endre kappa def til arctan ? (se:'destabilizing..') 
//           

// GSL Manual 26.6
//  http://stackoverflow.com/questions/27913858/how-to-consult-gsl-ode-with-many-parameters-and-harmonic-functions
// http://www.tutorialspoint.com/matlab/matlab_differential.htm
// Løser differentialene vha octave-symbolic se ~/xoct/minsky/differentiate.
// Transformerer octave-løsningene til c++ (e)^ -> exp() og x^y -> pow(x,y)
// se testdiff.cpp i ~/xlin/cpp/scratch
// rk8 bruker ikke jac fjernet f.o.m versjon 014
//
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <fstream>
#include <iomanip> // setiosflags(ios::fixed)
#include <limits>
#include <map>
#include <string>
#include <sstream>
#include <vector>

#include <boost/array.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/foreach.hpp>
#include <unistd.h>

using namespace std;

// Warn about use of deprecated functions.
#define GNUPLOT_DEPRECATE_WARN
#include "gnuplot-iostream.h"

#ifndef M_PI
#       define M_PI 3.14159265358979323846
#endif

inline void mysleep(unsigned millis){
  ::usleep(millis * 1000);
}

//http://stackoverflow.com/questions/7248627/setting-width-in-c-output-stream:
class formatted_output
{
private:
  int width, precision;
  std::ostream& stream_obj;
public:
  formatted_output(std::ostream& obj, int w, int p): width(w), precision(p),stream_obj(obj) {}
  template<typename T>
  formatted_output& operator<<(const T& output)  {
    stream_obj << std::fixed << std::setw(width) << std::setprecision(precision) << output;
    return *this;
  }
  formatted_output& operator<<(std::ostream& (*func)(std::ostream&))  {
    func(stream_obj);
    return *this;
  }
};

struct param_type {
  int eqset;
  double alfa;
  double beta;
  double delta;
  double nu;
  double r;
  double k1;
  double k2; 
  double k3;
  double k4;
  double k5;
  double k6;
  double k7;
  double k8;
} ;

double kappa(double k1, double k2, double k3, double profit){
  double response;
  response = k1 + exp(k2) * exp( k3 * profit);
  return response;
}  


double phi(double k4, double empshare){
  double response;
  response = (pow(k4,3)/(1.0-pow(k4,2)))/pow((1.0 - empshare),2) - k4/(1.0-pow(k4,2));
  return response;
}


double psi(double k5, double k6, double k7, double k8, double growthrate){
  double response;
  response = k5 + k6* exp(k7)*exp(k8 * growthrate);
  return response;
}

int func (double t, const double y[], double f[],
           void *params)
     {
       struct param_type *my_params_pointer = (param_type*) params;
        int eqset     = my_params_pointer->eqset;
	double alfa   = my_params_pointer->alfa;
	double beta   = my_params_pointer->beta;
	double delta  = my_params_pointer->delta;
	double nu     = my_params_pointer->nu;
	double r      = my_params_pointer->r;
	double k1     = my_params_pointer->k1;
	double k2     = my_params_pointer->k2;
	double k3     = my_params_pointer->k3;
	double k4     = my_params_pointer->k4;
	double k5     = my_params_pointer->k5;
	double k6     = my_params_pointer->k6;
	double k7     = my_params_pointer->k7;
	double k8     = my_params_pointer->k8;

        if (eqset == 74) {
	  f[0] =  y[0] * ( phi(k4, y[1]) - alfa) ;
	  f[1] =  y[1] * ( kappa(k1, k2, k3, 1 - y[0] - r * y[2]) /nu  - alfa - beta -delta);
	  f[2] =  y[2] * ( r -  kappa(k1, k2, k3, 1 - y[0] - r* y[2])/nu + delta) + kappa(k1, k2, k3,1 - y[0] - r * y[2]) - (1 - y[0]) + y[3];
	  f[3] =  y[3] * ( psi(k5, k6, k7, k8, kappa(k1, k2, k3, 1 - y[0] -r*y[2]) / nu  - delta) + kappa(k1, k2, k3, 1 - y[0] - r * y[2]) / nu - delta );
	}
	
	if (eqset == 77) {
	  f[0] =  y[0] * ( phi(k4, y[1]) - alfa) ;
	  f[1] =  y[1] * (( kappa(k1, k2, k3, 1 - y[0] - r/( y[2] * y[3] )) /nu ) - alfa - beta -delta);
	  f[2] =  y[2] * ( psi(k5, k6, k7, k8, kappa(k1, k2, k3, 1 - y[0] -r/(y[2] * y[3] )) / nu - delta) - r) - ( y[2] * y[2] * ( y[3] * kappa(k1, k2, k3,1 - y[0] - r/(y[2] * y[3])) - y[3] * (1 - y[0]) + 1 ));
	  f[3] =  y[3] * ( -psi(k5, k6, k7, k8, ( kappa(k1, k2, k3, 1 - y[0] -r/(y[2] * y[3])) / nu ) - delta) + (kappa(k1, k2, k3, 1 - y[0] - r/( y[2] * y[3] )) / nu) - delta );
	}

        return GSL_SUCCESS;
     }

     // fjernet - finnes i minsky.jac.cc = minsky.013.cc    
     /*
     int jac (double t, const double y[], double *dfdy, 
          double dfdt[], void *params)
     {
      return GSL_SUCCESS;
     }
     */


int main(int argc, char **argv) {

       if (argc < 2) {
         std::cout << "Usage: minsky eqset (74|77) [omega lambda d|upsilon  p|x num_steps])" << endl;
         exit(0);
       }

       int eqset = 74;	 
       if (argc > 1) eqset   = atoi(argv[1]);
       cout << eqset << endl; // CONTROL
       if (eqset != 74 && eqset != 77) {
         std::cout << "eqset has to be either 74 or 77" << endl;
         exit(0);
       }
       double y[4] = { 0.0, 0.0, 0.0, 0.0 } ; 
       if (eqset == 74){
	 //y[4] = { 0.95, 0.9, 0.0, 0.001 } ;   // omega, lambda, d, p
	 y[0] = 0.95;
	 y[1] = 0.9;
	 y[2] = 0.0;
	 y[3] = 0.001;
       }
       if (eqset == 77) {
	 //y[4] = { 0.95, 0.9, 0.12, 100 } ;   // omega, lambda, upsilon, x
	 y[0] = 0.95;
	 y[1] = 0.9;
	 y[2] = 0.12;
	 y[3] = 100.0;
	 
       }
       int num_steps = 100;
         
       if (argc > 2) y[0] = atof(argv[2]);
       if (argc > 3) y[1] = atof(argv[3]);
       if (argc > 4) y[2] = atof(argv[4]);
       if (argc > 5) y[3] = atof(argv[5]);
       if (argc > 6) num_steps = atoi(argv[6]);

       Gnuplot gp1;
       Gnuplot gp2;
       Gnuplot gp3;
       Gnuplot gp4;      

       std::ostringstream oss;
       oss << "set title \"";
       for (int j = 0; j < argc; j++)
	 { oss << argv[j] << ' '; }
       oss << "\\n\"" << std::endl ;
       std::cout << oss.str() << std::endl;

       std::ofstream funkout("minsky.funk.out");
       std::ofstream varout("minsky.var.out");       
       formatted_output fout(funkout, 16, 6);
       formatted_output vout(varout, 12, 6);              

       std::vector<std::pair<double, double> > xy_pts_A; // omega
       std::vector<std::pair<double, double> > xy_pts_B; // lambda
       std::vector<std::pair<double, double> > xy_pts_C; // d         // upsilon
       std::vector<std::pair<double, double> > xy_pts_D; // p         // x

       std::vector<std::pair<double, double> > xy_pts_E; // p (ponzi)      
       std::vector<std::pair<double, double> > xy_pts_F; // d (debt)
       std::vector<std::pair<double, double> > xy_pts_G; // g (growth)
       std::vector<std::pair<double, double> > xy_pts_Y; // Y

       double alfa   = 0.025;
       double beta   = 0.02;
       double delta  = 0.01;
       double nu     = 3.0;
       double r      = 0.03;
       double k1     = -0.0065;
       double k2     = -5 ;
       double k3     =  20;
       double k4     = 0.04;
       double k5     = -0.25;
       double k6     = 0.25;
       double k7     = -0.36;
       double k8     = 12.0;
       double growth = 0.0;
       double Y      = 0.0;

       double zphi;
       double zkappa;
       double zpsi;
       
       if (eqset == 74) {
         growth = kappa(k1, k2, k3, 1 - y[0] - r * y[2]) / nu - delta;
       }
       else if (eqset == 77) {
	 growth = kappa(k1, k2, k3, 1 - y[0] - r/( y[2] * y[3] )) / nu - delta;
       }
       Y = 100.0 * ( 1.0 + growth);
       
       struct param_type my_params = {eqset, alfa, beta, delta, nu, r, k1, k2, k3, k4, k5, k6, k7, k8 };

       // gsl_odeiv2_system sys = {func, jac, 4, &my_params};
       // http://www.physics.buffalo.edu/phy411-506/tools/gsl/ode/index.html
       // rk8pd bruker ikke jac settes til NULL:
       // gsl_odeiv2_system sys = {func, jac, 4, &my_params};
       gsl_odeiv2_system sys = {func, NULL, 4, &my_params};
       gsl_odeiv2_driver * d = 
       gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd,
					1e-6, 1e-6, 0.0);
       int i = 0;

       double t = 0.0, t1 = 100.0;

       xy_pts_A.push_back(std::make_pair(t, y[0]));             //omega
       xy_pts_B.push_back(std::make_pair(t, y[1]));		//lambda	    
       xy_pts_C.push_back(std::make_pair(t, y[2]));             //d|upsilon
       xy_pts_D.push_back(std::make_pair(t, y[3]));		//p|x
       
       if (eqset == 77) {
	 xy_pts_E.push_back(std::make_pair(t, 1 / y[3]));         //ponzi
	 xy_pts_F.push_back(std::make_pair(t, 1 / y[3] / y[2]));  //debt
       }
       xy_pts_G.push_back(std::make_pair(t, growth));           //growth
       xy_pts_Y.push_back(std::make_pair(t, Y));                //Y

       zphi   = phi(k4, y[1]);
       if (eqset == 74) {
	 zkappa = kappa(k1, k2, k3, 1 - y[0] - r * y[2]);
	 zpsi   = psi(k5, k6, k7, k8, kappa(k1, k2, k3, 1 - y[0] -r * y[2]) / nu - delta);
       }
       if (eqset == 77) {
	 zkappa = kappa(k1, k2, k3, 1 - y[0] - r/( y[2] * y[3] ));
	 zpsi   = psi(k5, k6, k7, k8, kappa(k1, k2, k3, 1 - y[0] -r/(y[2] * y[3] )) / nu - delta);
       }
       fout << oss.str() << std::endl;
       fout << "t" << "phi" <<"kappa" << "psi" << std::endl;
       fout << t << zphi << zkappa << zpsi << std::endl;       

       vout << oss.str() << std::endl;
       if (eqset == 74)  {
	 vout << "t" << "omega" <<"lambda" << "d" << "p" << "growth" << "Y"  << std::endl;
	 vout << t  << y[0] << y[1] << y[2] << y[3] << growth << Y <<std::endl;       
       }
       if (eqset == 77) {
 	 vout << "t" << "omega" <<"lambda" << "upsilon" << "debt" << "ponzi" << "x" << "growth" << "Y"  << std::endl;
	 vout << t  << y[0] << y[1] << y[2] << y[3] <<  1/y[3]/y[2] <<1/y[3] << growth << Y <<std::endl;       }
       if (eqset == 74) {
	 printf ("%6s %10s %10s %10s %10s %10s %10s %10s %10s\n", "t","omega","lambda","d","p","","","growth","Y");
	 printf ("%6.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n", t, y[0], y[1], y[2], y[3], 0.0, 0.0, growth, Y);
       }
       if (eqset == 77) {
         printf ("%6s %10s %10s %10s %10s %10s %10s %10s %10s\n", "t","omega","lambda","upsilon","x","ponzi","debt","growth","Y");
         printf ("%6.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n", t, y[0], y[1], y[2], y[3], 1/y[3]/y[2], 1/y[3],growth, Y);
       }

       
       for (i = 1; i <= num_steps; i++)
         {
           double ti = i * t1 / 100.0;
           int status = gsl_odeiv2_driver_apply (d, &t, ti, y);
           if (status != GSL_SUCCESS)
	     {
	       printf ("error, return value=%d\n", status);
	       break;
	     }

           //Funksjonsverdiene:
	   zphi   = phi(k4, y[1]);
	   if (eqset == 74) {
	     zkappa = kappa(k1, k2, k3, 1 - y[0] - r * y[2]);
	     zpsi   = psi(k5, k6, k7, k8, kappa(k1, k2, k3, 1 - y[0] -r *y[2]) / nu - delta);
	     growth = kappa(k1, k2, k3, 1 - y[0] - r* y[2]) / nu - delta;
	   }
	   if (eqset == 77) {
	     zkappa = kappa(k1, k2, k3, 1 - y[0] - r/( y[2] * y[3] ));
	     zpsi   = psi(k5, k6, k7, k8, kappa(k1, k2, k3, 1 - y[0] -r/(y[2] * y[3] )) / nu - delta);
	     growth = kappa(k1, k2, k3, 1 - y[0] - r/( y[2] * y[3] )) / nu - delta;
	   }
	   Y = Y * ( 1 + growth);

	   fout << t << zphi << zkappa << zpsi << std::endl;
	   vout << t << y[0] << y[1] << y[2] << y[3] << growth << Y <<std::endl;       	   
	   if (eqset == 74) {
	     printf ("%6.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n", t, y[0], y[1], y[2], y[3], 0.0, 0.0, growth, Y);
	   }
	   if (eqset == 77) {
	     printf ("%6.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n", t, y[0], y[1], y[2], y[3], 1/y[3]/y[2], 1/y[3], growth, Y);
	   }
 	   xy_pts_A.push_back(std::make_pair(t, y[0]));
	   xy_pts_B.push_back(std::make_pair(t, y[1]));			    
	   xy_pts_C.push_back(std::make_pair(t, y[2]));
	   xy_pts_D.push_back(std::make_pair(t, y[3]));			    
	   if (eqset == 77) {
	     xy_pts_E.push_back(std::make_pair(t, 1 / y[3]));         //ponzi
	     xy_pts_F.push_back(std::make_pair(t, 1 / y[3] / y[2]));  //debt
	   }
	   xy_pts_G.push_back(std::make_pair(t, growth)); //growth
           xy_pts_Y.push_back(std::make_pair(t, Y)); //Y
	 }
     
       gsl_odeiv2_driver_free (d);

       // GNPLOTS:
       if (eqset == 74){
	 // Omega & Lambda:
	 gp1 << "set terminal x11 1 persist size 640,450 position 50,50\n";
	 gp1 << "set grid\n";
	 gp1 << oss.str() ; // Title string
	 gp1 << "set ylabel 'omega lambda'\n";
	 gp1 << "plot '-' with lines title 'omega', '-' with lines title 'lambda'\n";
	 gp1.send1d(xy_pts_A);
	 gp1.send1d(xy_pts_B);
	 // p & d:
	 gp2 << "set terminal x11 1 persist size 640,450 position 50,536\n";
	 gp2 << "set grid\n";
	 gp2 << oss.str() ; // Title string
	 gp2 << "set ylabel 'p\n";
	 gp2 << "set y2label 'd'\n";
	 gp2 << "set y2tics\n";
	 //gp2 << "plot '-' with lines title 'd\n";
 	 gp2 << "plot '-' with lines title 'p', '-' with lines  title 'd' axes x1y2\n";
	 gp2.send1d(xy_pts_D);
	 gp2.send1d(xy_pts_C);
	 
	 /*	 // p:
	 gp3 << "set terminal x11 1 persist size 640,450 position 690,536\n";
	 gp3 << "set grid\n";
	 gp3 << oss.str() ; // Title string
	 gp3 << "set ylabel 'p\n";
	 gp3 << "plot '-' with lines title 'p\n";
	 gp3.send1d(xy_pts_D);
	 */
       }

       if (eqset == 77) {
	 // Omega & Lambda & Upsilon & X:
	 gp1 << "set terminal x11 1 persist size 640,450 position 50,50\n";
	 gp1 << "set grid\n";
	 gp1 << oss.str() ; // Title string
	 gp1 << "set ylabel 'omega lambda upsilon'\n";
	 gp1 << "set y2label 'x'\n";
	 gp1 << "set y2tics\n";
	 gp1 << "plot '-' with lines title 'omega', '-' with lines title 'lambda', '-' with lines title 'upsilon','-' with lines  title 'x' axes x1y2\n";
	 gp1.send1d(xy_pts_A);
	 gp1.send1d(xy_pts_B);
	 gp1.send1d(xy_pts_C);
	 gp1.send1d(xy_pts_D);
	 // Ponzi & Debt: 
	 gp2 << "set terminal x11 2 persist size 640,450 position 50,536\n";
	 gp2 << "set grid\n";
	 gp2 << oss.str() ; // Title string
	 gp2 << "set ylabel 'ponzi'\n";
	 gp2 << "set y2label 'debt'\n";
	 gp2 << "set y2tics\n";
	 gp2 << "plot '-' with lines title 'ponzi', '-' with lines title 'debt' axes x1y2\n";
	 gp2.send1d(xy_pts_E);
	 gp2.send1d(xy_pts_F);
       }		    
       //  p/ponzi & growth:
       gp3 << "set terminal x11 2 persist size 640,450 position 690,50\n";
       gp3 << "set grid\n";
       gp3 << oss.str() ; // Title string
       gp3 << "set ylabel 'p/ponzi'\n";
       gp3 << "set y2label 'growthi'\n";
       gp3 << "set y2tics\n";
       gp3 << "plot '-' with lines title 'p/ponzi', '-' with lines title 'growth' axes x1y2\n";
       if (eqset == 74) gp3.send1d(xy_pts_D);
       if (eqset == 77) gp3.send1d(xy_pts_E);              
       gp3.send1d(xy_pts_G);        

       //  Y & d/debt:
       gp4 << "set terminal x11 2 persist size 640,450 position 690,536\n";
       gp4 << "set grid\n";
       gp4 << oss.str() ; // Title string
       gp4 << "set ylabel 'Y'\n";
       gp4 << "set y2label 'd/debt'\n";
       gp4 << "set y2tics\n";
       gp4 << "plot '-' with lines title 'Y', '-' with lines title 'd/debt' axes x1y2\n";
       gp4.send1d(xy_pts_Y);       
       if (eqset == 74) gp4.send1d(xy_pts_C);
       if (eqset == 77) gp4.send1d(xy_pts_F);              

       return 0;    
     }

// EOF

Tags: