function model = gro2model()

symbols = struct;
symbols.states = {'K','Z'};
symbols.controls = {'C','I','Mu'};
symbols.expectations = {'E'};
symbols.shocks = {'Epsilon'};
symbols.parameters = {'a','tau','delta','beta','rho','alpha'};

calibration = struct;
calibration.states = [ 1.  0.];
calibration.controls = [ 0.22117193  0.          0.        ];
calibration.expectations = [ 4.75933623];
calibration.parameters = [ 0.24077193  1.          0.0196      0.95        0.9         0.3       ];
calibration.sigma = [[ 0.]];

functions = struct;
functions.transition = @transition;
functions.arbitrage = @arbitrage;
functions.expectation = @expectation;

functions.arbitrage_lb = @arbitrage_lb;
functions.arbitrage_ub = @arbitrage_ub;

model = struct;
model.symbols = symbols;
model.functions = functions;
model.calibration = calibration;
model.infos.incidence_matrices.arbitrage_lb = {[0 0 ; 0 0 ; 0 0] };
model.infos.incidence_matrices.arbitrage_ub = {[0 0 ; 0 0 ; 0 0] };
model.infos.incidence_matrices.transition = {[1 0 ; 0 1] [0 1 0 ; 0 0 0] [0 ; 1] };
model.infos.incidence_matrices.arbitrage = {[1 1 ; 0 0 ; 0 0] [1 1 0 ; 1 0 1 ; 0 0 1] [0 ; 0 ; 1] };
model.infos.incidence_matrices.expectation = {[0 0] [0 0 0] [0] [1 1] [1 0 1] };

end

function [val, val_states_p, val_controls_p, val_shocks] = transition(states_p, controls_p, shocks, p, output)

    if nargin <= 4
        output = zeros(4,1);
        for i = 1:nargout
            output(i) = 1;
        end
    end

    n = size(states_p,1);

    val = zeros( n, 2 );
    val(:,1) = controls_p(:,2) + states_p(:,1).*(-p(3) + 1);
    val(:,2) = p(5).*states_p(:,2) + shocks(:,1);

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: states_p

    if output(2)
        val_states_p = zeros( n,2,2 );
        val_states_p(:,1,1) = -p(3) + 1;
        val_states_p(:,2,2) = p(5);
    else
        val_states_p = [];
    end;

    % Derivatives w.r.t: controls_p

    if output(3)
        val_controls_p = zeros( n,2,3 );
        val_controls_p(:,1,2) = 1;
    else
        val_controls_p = [];
    end;

    % Derivatives w.r.t: shocks

    if output(4)
        val_shocks = zeros( n,2,1 );
        val_shocks(:,2,1) = 1;
    else
        val_shocks = [];
    end;

end

function [val, val_states, val_controls, val_expectations] = arbitrage(states, controls, expectations, p, output)

    if nargin <= 4
        output = zeros(4,1);
        for i = 1:nargout
            output(i) = 1;
        end
    end

    n = size(states,1);

    val = zeros( n, 3 );
    val(:,1) = -p(1).*states(:,1).^p(6).*exp(states(:,2)) + controls(:,1) + controls(:,2);
    val(:,2) = controls(:,3) + controls(:,1).^(-p(2));
    val(:,3) = p(4).*expectations(:,1) + controls(:,3);

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: states

    if output(2)
        val_states = zeros( n,3,2 );
        val_states(:,1,1) = -p(1).*p(6).*states(:,1).^(p(6) - 1).*exp(states(:,2));
        val_states(:,1,2) = -p(1).*states(:,1).^p(6).*exp(states(:,2));
    else
        val_states = [];
    end;

    % Derivatives w.r.t: controls

    if output(3)
        val_controls = zeros( n,3,3 );
        val_controls(:,1,1) = 1;
        val_controls(:,1,2) = 1;
        val_controls(:,2,1) = -p(2).*controls(:,1).^(-p(2) - 1);
        val_controls(:,2,3) = 1;
        val_controls(:,3,3) = 1;
    else
        val_controls = [];
    end;

    % Derivatives w.r.t: expectations

    if output(4)
        val_expectations = zeros( n,3,1 );
        val_expectations(:,3,1) = p(4);
    else
        val_expectations = [];
    end;

end

function [val, val_states, val_controls, val_shocks_f, val_states_f, val_controls_f] = expectation(states, controls, shocks_f, states_f, controls_f, p, output)

    if nargin <= 6
        output = zeros(6,1);
        for i = 1:nargout
            output(i) = 1;
        end
    end

    n = size(states,1);

    val = zeros( n, 1 );
    val(:,1) = p(1).*p(6).*controls_f(:,1).^(-p(2)).*states_f(:,1).^(p(6) - 1).*exp(states_f(:,2)) - controls_f(:,3).*(-p(3) + 1);

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: states

    if output(2)
        val_states = zeros( n,1,2 );
    else
        val_states = [];
    end;

    % Derivatives w.r.t: controls

    if output(3)
        val_controls = zeros( n,1,3 );
    else
        val_controls = [];
    end;

    % Derivatives w.r.t: shocks_f

    if output(4)
        val_shocks_f = zeros( n,1,1 );
    else
        val_shocks_f = [];
    end;

    % Derivatives w.r.t: states_f

    if output(5)
        val_states_f = zeros( n,1,2 );
        val_states_f(:,1,1) = p(1).*p(6).*controls_f(:,1).^(-p(2)).*states_f(:,1).^(p(6) - 2).*(p(6) - 1).*exp(states_f(:,2));
        val_states_f(:,1,2) = p(1).*p(6).*controls_f(:,1).^(-p(2)).*states_f(:,1).^(p(6) - 1).*exp(states_f(:,2));
    else
        val_states_f = [];
    end;

    % Derivatives w.r.t: controls_f

    if output(6)
        val_controls_f = zeros( n,1,3 );
        val_controls_f(:,1,1) = -p(1).*p(6).*p(2).*controls_f(:,1).^(-p(2) - 1).*states_f(:,1).^(p(6) - 1).*exp(states_f(:,2));
        val_controls_f(:,1,3) = p(3) - 1;
    else
        val_controls_f = [];
    end;

end

function [val, val_s] = arbitrage_lb(s, p, output)

    if nargin <= 2
        output = zeros(2,1);
        for i = 1:nargout
            output(i) = 1;
        end
    end

    n = size(s,1);

    val = -inf( n , 3 );
    val(:,2) = 0;

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: s

    if output(2)
        val_s = zeros( n,3,2 );
    else
        val_s = [];
    end;

end

function [val, val_s] = arbitrage_ub(s, p, output)

    if nargin <= 2
        output = zeros(2,1);
        for i = 1:nargout
            output(i) = 1;
        end
    end

    n = size(s,1);

    val = inf( n , 3 );

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: s

    if output(2)
        val_s = zeros( n,3,2 );
    else
        val_s = [];
    end;

end