# User:Jvgama

% Obtains the optimal intervals N1, ..., NJ. % M(n) > 0, K > 0. % From Syst04Pc.

% Permanent shock case (eta = 0)

% January 26, 2015

function y = SystPerm01c(x)

global gamma; global M0n; global Nf; global P0; global P; global r2; global rho; global Y;

L = length(P); Lx = length(x); method = 'linear';  % Interpolation method. incr = 0.1; lambda = 1/(P0*Y);

TrTj = cumsum(x);  % Transfer times, Tj. PTj = interp1(0:(L-1),P,TrTj,method);  % Price at time Tj. PTjp = interp1(0:(L-1),P,TrTj+incr,method);  % Price at time Tj+incr. InfTj = (log(PTjp)-log(PTj))/incr;  % Inflation at Tj.

% Case eta > 0 %rt = r1 + (r2-r1)*exp(-eta*TrTj);  % 1xNNj %Rt = r1*TrTj - (r2-r1)/eta*exp(-eta*TrTj) + (r2-r1)/eta;  % 1xNNj

% Case eta = 0 => r(t) = r2 and R(t) = r2*t Rt = r2*TrTj;  % 1xNNj

miuL = exp(-Rt(1))*lambda; KnL = M0n - (1-exp(-rho*TrTj(1)))/rho/miuL;  % Kn > 0.

% Matrix y % 1st entry (A, B, C): column 1 % 2nd entry (D, E, F): columns 2 to NNj-1 % 3rd entry (G, H, IA): column NNj

A = exp(-Rt(1)+rho*TrTj(1))*gamma*PTj(1)*(r2-InfTj(1))/P0; B = r2*(1-exp(-rho*x(2)))/rho; C = -exp(-Rt(1)+rho*TrTj(1))*r2*KnL*lambda;

D = -r2*x(2:(Lx-1)); E = gamma*PTj(2:(Lx-1)).*exp(-Rt(2:(Lx-1))+rho*TrTj(2:(Lx-1))).*(r2-InfTj(2:(Lx-1)))/P0; F = r2*(1-exp(-rho*x(3:Lx)))/rho;

G = -r2*x(Lx); H = gamma*PTj(Lx)*exp(-Rt(Lx)+rho*TrTj(Lx))*(r2-InfTj(Lx))/P0; IA = r2*(1-exp(-rho*Nf))/rho;

y = [A + B + C,...

```    D + E + F,...
G + H + IA];
```