The Goodwin Model is essentially a 2 equation predator-prey system.

The first equation is for u, the wage share. For nominal wages w,
the wage bill for nominal output Y = wL/Y = wL/paL = w/pa .

Writing u’ for du/dt then the first equation is

(1) u’/u = w’/w – a’/a -p’/p

w’/w = Ph(λ) is the Phillips Curve relationship for wage growth based on employment level.

a’/a = α , the productivity growth rate is also taken as constant α

p’/p = g_p is the assumed piecewise constant inflation rate.

(1) then becomes (1a) u’/u = Ph(λ) – α – g_p is the final form of the wage share growth rate equation.

From the relation Q = aL we can write
Q’/Q = α + L’/L = α + λ’/λ + n’/n

Since L = λn where n is the population and λ is the employment percentage. n’/n = β is assumed piecewise constant so Q’/Q = λ’/λ + α + β.

From Q = k/σ, Q’/Q = k’/k

Also, k’ = b(1-u)Q – δk so that Q’/Q can be written

b(1-u)Q/σQ – δ = b(1-u)/σ – δ

δ is a constant capital depreciation rate.

Finally, regrouping

(2) λ’/λ = b(1-u)/σ -(α + β + δ) is the employment percentage growth rate.

Scilab Code

z0=[.97;.95;2;1;25;8.33 ]; //[wage share;%employment;wage;price level; real GDP billions]
t=1913:.2:2100;
clf()

function dx =minsky1(t,x)

alfa=.025//+.05*(x(1)-.9);  // alfa is productivity growth rate
g_p=0     //g_p is exogenous inflation rate. Source is measuringworth.com
bta=.0123//+.05*(x(1)-.9)//  bta is population growth rate
dlta=.01 //  dlta is capital depreciation rate

if t>1940 then
g_p=.040
end

ex=0

//if t>1929 then
//ex=0
//end
//
//if t>1939 then
//ex=1
//end
//
//ex=0

v=.33*x(4)*x(2)^ex  // v is constant accelerator multiplied by price level
//v=.33

if t > 1979 then
g_p=.00
end

if t>1984 then
g_p=.04
end

T0 = 2112//97629595999999
B=.1
C=5

if t>T0 then  //Return of Gold Standard
g_p=-.0123
end

//if t>T0 then
//g_p=-bta//*x(2)^(ex)
//end

//g_p=0

b=1//*(1-g_p)// b is investment fraction K'/K = (1-u)b/v-dlta  etc.

function GE=genexp(x,xv,yv,s,m)  // Phillips curve function
GE=(yv-m)*exp(s/(yv-m)*(x-xv))+m
endfunction

Ph=genexp(x(2),.9,0,.8,-.05)
n1=100

dx(1)=x(1)*(1-2*x(1)^n1/(n1+2))^.5*(Ph-alfa-g_p)
//dx(1)=x(1)*(Ph-alfa-g_p) // x(1) = worker's share of output

//dx(2)=x(2)*b*((1-x(1))/v-alfa-bta-dlta-g_p) // x(2) = %employment
dx(2)=x(2)*(1-2*x(2)^n1/(n1+2))^.5*b*((1-x(1))/v-alfa-bta-dlta)
dx(3)=x(3)*(Ph+B*(-exp(-(t-T0)^2/.15)-exp(-(t-(T0+C))^2/.15)-exp(-(t-(T0+2*C))^2/.15)))//%employment

dx(4)=x(4)*g_p+1*exp(-(t-1933)^2/.17)-x(4)*B*exp(-(t-T0)^2/.2)-x(4)*B*exp(-(t-(T0+C))^2/.2)-x(4)*B*exp(-(t-(T0+2*C))^2/.2)// price level 

//dx(5)=(alfa+bta*x(2))*x(5)  // Real Output
dx(5)=(b*(1-x(1))/v-dlta)*x(5)  // Real Output

dx(6)=x(6)*(b*(1-x(1))/v-dlta)  //Real capital

endfunction

z=ode(z0,1913,t,minsky1);

plot2d(t',[z(1,:);z(2,:);z(4,:)^(-1)]',style=[1 2 5]);
legends(['wage share';'%employment';'1/p'],[1 2 5],"ll")
halt
clf();
plot2d4(z(1,:)',z(2,:)');
legends(' %employment vs wage share ',style=1,"ll");
halt
clf();

plot2d(t',[z(3,:);z(3,:)./z(4,:);z(4,:)]',logflag="nl",style=[1 2 3 4 ]);
legends(['wages';'real wages';'price level'],[1 2 3 ],"ul");
halt
clf
y=z(5,:);
py=z(4,:).*y;

plot2d(t',[y;py;py./y;z(6,:)]',[1 2 5 3],logflag="nl")
legends(['Q = real output';'Y = nominal output';'Y/Q ratio';'K = real Capital'],[1 2 5 3],"ul")