ShapeAttack/main.m at master · kellinpelrine/ShapeAttack · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
clear
close all
rng(12378)

tic
options = optimset('Display','iter');

gamma       = 1;                %assuming log utility
epsilon     = 2;                %elasticity of labor supply
alpha       = .4;               %capital's share of income
rho         = .95;
sigmaerror  =.01;
mu          = 0;                %long-run growth rate
%Steady State Values
IYssratio   = .25;
KYssratio   = 14;
nss         = 1/3;
kss         = 14^(.06) / 3;
%Calibrated Parameter Values
delta   = 1/56;                  %depreciation rate
beta    = .964;                  %discount rate
b       = 7.2;                  %disutility of labor multiple

fun = @(x) -abs(model(gamma,epsilon,alpha,rho,sigmaerror,mu,IYssratio,KYssratio,nss,kss,delta,x,b) - .745);
%[test, argtest] = fminbnd(fun,.8,.999,options)
%test = model(gamma,epsilon,alpha,rho,sigmaerror,mu,IYssratio,KYssratio,nss,kss,delta,beta,b)

balllow = .5;
ballhigh = .99;
stepsize = .02;
mutationprob = .5;
popsize = 12;
beta0 = .964;
generations = 10;

%bernp = rand(popsize,1) <= mutationprob;
%x = beta0*ones(popsize,1) + bernp .* (  rand(popsize,1)*( (ballhigh - balllow)*stepsize ) - balllow*stepsize  ) ;
% %x = beta0*ones(popsize,1) + (  rand(popsize,1)*( (ballhigh - balllow)*stepsize ) - balllow*stepsize  ) ;

%x = linspace(balllow + .05,ballhigh - .05,popsize)' + (  rand(popsize,1)*( (ballhigh + balllow)*stepsize ) - balllow*stepsize  ) ;
x = linspace(.8,.95,popsize)' + (  rand(popsize,1)*( (ballhigh - balllow)*stepsize ) - balllow*stepsize  ) ;
%x = linspace(balllow + .05,ballhigh - .05,popsize)' + (  randn(popsize,1) * sqrt(stepsize) ) ;
F = zeros(popsize,generations);
children = zeros(popsize - 1,1);
parents = zeros(popsize - 1,1);
xsave = zeros(popsize,generations); %for analyzing output
figure(1)
hold on
for g = 1:generations
    x(x < balllow) = balllow;
    x(x > ballhigh) = ballhigh;
    xsave(:,g) = x;
    plot(xsave(:,g))
    parfor p = 1:popsize
        F(p,g) = -fun(x(p));
    end
    [Fmax,Fargmax] = max(F(:,g));
    %Fmax
    g
    %x(Fargmax)
    probs = F(:,g) ./ sum(F(:,g));
    parents = randsrc(popsize - 1,2,[linspace(1,popsize,popsize) ; probs']);
    children(1:popsize - 1) = ( x(parents(:,1)).*F(parents(:,1),g) + x(parents(:,2)).*F(parents(:,2),g) ) ./ ( F(parents(:,1),g) + F(parents(:,2),g) );
    mutate = rand(popsize - 1,1) <= mutationprob;
    %children = children + mutate .* (  rand(popsize - 1,1)*( (ballhigh + balllow)*stepsize ) - balllow*stepsize  );
    children = children + mutate .* (  randn(popsize - 1,1) * sqrt(stepsize) ) ;
    x(1) = x(Fargmax);
    x(2:popsize) = children;
end

hold off
plot(xsave')
elitefitness = max(F,[],1);
figure(2)
plot(elitefitness)
toc