function model = cs1model()

symbols = struct;
symbols.states = {'X'};
symbols.controls = {'C'};
symbols.expectations = {'E'};
symbols.shocks = {'Y'};
symbols.parameters = {'r','delta','rho'};

calibration = struct;
calibration.states = [ 100.];
calibration.controls = [ 100.];
calibration.expectations = [ 0.0001];
calibration.parameters = [ 0.05  0.1   2.  ];
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] };
model.infos.incidence_matrices.arbitrage_ub = {[1] };
model.infos.incidence_matrices.transition = {[1] [1] [1] };
model.infos.incidence_matrices.arbitrage = {[0] [1] [1] };
model.infos.incidence_matrices.expectation = {[0] [0] [0] [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, 1 );
    val(:,1) = shocks(:,1) + (p(1) + 1).*(-controls_p(:,1) + states_p(:,1));

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: states_p

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

    % Derivatives w.r.t: controls_p

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

    % Derivatives w.r.t: shocks

    if output(4)
        val_shocks = zeros( n,1,1 );
        val_shocks(:,1,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, 1 );
    val(:,1) = expectations(:,1).*(p(1) + 1)./(p(2) + 1) - controls(:,1).^(-p(3));

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: states

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

    % Derivatives w.r.t: controls

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

    % Derivatives w.r.t: expectations

    if output(4)
        val_expectations = zeros( n,1,1 );
        val_expectations(:,1,1) = (p(1) + 1)./(p(2) + 1);
    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) = controls_f(:,1).^(-p(3));

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: states

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

    % Derivatives w.r.t: controls

    if output(3)
        val_controls = zeros( n,1,1 );
    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,1 );
    else
        val_states_f = [];
    end;

    % Derivatives w.r.t: controls_f

    if output(6)
        val_controls_f = zeros( n,1,1 );
        val_controls_f(:,1,1) = -p(3).*controls_f(:,1).^(-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 , 1 );

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: s

    if output(2)
        val_s = zeros( n,1,1 );
    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 , 1 );
    val(:,1) = s(:,1);

    if nargout <= 1
        return
    end

    % Derivatives w.r.t: s

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

end