%--------------------------------------------------------------------------
% PARAMETER CONTROLING MODEL CREATION
%--------------------------------------------------------------------------
n = 500; % number of endogenous variables
m = 50; % minimal number of state variables (other may appear)
k = 10; % number of fundamental shocks (sunspot shocks may appear)
min = -.5; % minimal predefined eigenvalue of state transition matrix
max = 0.9; % maximal predefined eigenvalue of state transition matrix
xi = 1; % growth restriction;
do_reduction = 1; % do model reduction
method = 1; % method of calculating null spaces
disp(' ');
disp('Model:');
disp([' ' num2str(n) ' - number of endogenous variables']);
disp([' ' num2str(m) ' - minimal number of state variables (other may appear)']);
disp([' ' num2str(k) ' - number of fundamental shocks (sunspot shocks may appear)']);
disp(' ');
disp([' ' num2str(min) ' - minimal predefined eigenvalue of state transition matrix']);
disp([' ' num2str(max) ' - maximal predefined eigenvalue of state transition matrix']);
disp([' ' num2str(xi) ' - growth restriction']);
disp(' ');
disp([' ' num2str(do_reduction) ' - do model reduction']);
disp([' ' num2str(method) ' - method of calculating null spaces']);
%--------------------------------------------------------------------------
% MODEL CREATION
%--------------------------------------------------------------------------
disp(' ');
disp('Preparing matrices');
B = [eye(k) zeros(k,n-k);zeros(n-k,n)];
C = full(sprand(n,n,1/n/10));
R = randn(n,m);
Q = randn(m,k);
S = randn(n,k);
a = (max-min)/(m-1);
U = orth(randn(m,m));
P = U*diag([min:a:max])*U';
temp = null(R');
A = (R'\(-(B+C)*R*P)')'+randn(n,size(temp,2))*temp';
V = -A*S;
W = -B*S+B*R*Q;
%--------------------------------------------------------------------------
% SOLVING THE MODEL
%--------------------------------------------------------------------------
disp(' ');
disp('Solving the model');
tic
[Pbar,Rbar,Sbar,Qbar,N1bar,N2bar,eig,info] = linear2ss(A,B,C,V,W,xi,do_reduction,method);
t1 = toc;
X1 = randn(size(N1bar,2),size(Sbar,2));
X2 = randn(size(N2bar,2),0);
Sbar1 = Sbar+N1bar*X1;
Sbar2 = N1bar*X2;
Qbar1 = Qbar+N2bar*X1;
Qbar2 = N2bar*X2;
n1 = norm(A*Rbar+(B+C)*Rbar*Pbar);
n2 = norm(A*Sbar1+V);
n3 = norm(A*Sbar2);
n4 = norm(B*(Rbar*Qbar1+Sbar1)+W);
n5 = norm(B*(Rbar*Qbar2+Sbar2));
n = norm([n1 n2 n3 n4 n5]);
disp(['Model solved in ' num2str(t1) ' seconds']);
disp(['Norm of residuals from equations describing solution to the model ' num2str(n) ]);
disp(' ');
disp('Solution:');
disp([' ' num2str(size(Pbar,1)) ' - number of state variables']);
disp([' ' num2str(size(N1bar,2)) ' - number of sunspot shocks']);
%--------------------------------------------------------------------------
% SETTING DIFFERENTIALS OF MODEL MATRICES WITH RESPECT TO
% PARAMETERS
%--------------------------------------------------------------------------
Sdif = zeros(size(Sbar)); Sdif(1,1) = 1;
Qdif = zeros(size(Qbar)); Qdif(1,2) = 1;
Adif = zeros(size(A)); Adif(1,1) = 1;
Bdif = zeros(size(B)); Bdif(1,1) = 1;
Cdif = zeros(size(C)); Cdif(1,1) = 1;
% creating differentials
Vdif = -A*Sdif-Adif*S;
Wdif = -Bdif*(Rbar*Qbar+Sbar)-B*(Rbar*Qdif+Sdif);
Ag = {A, Adif};
Bg = {B, Bdif};
Cg = {C, Cdif};
Vg = {V, Vdif};
Wg = {W, Wdif};
%--------------------------------------------------------------------------
% FIND DIFFERENTIALS
%--------------------------------------------------------------------------
disp(' ');
disp('Calculating differentials of solution matrices with respect to parameters');
tic
[Pdif,Rdif,S1dif,S2dif,Q1dif,Q2dif,INFO] = gradient(Ag,Bg,Cg,Vg,Wg,Pbar,Rbar,Sbar1,Sbar2,Qbar1,Qbar2,method,1);
t2=toc;
disp(['Differentials calculated in ' num2str(t2) ' seconds']);
n1 =norm(Ag{2}*Rbar+(Bg{2}+Cg{2})*Rbar*Pbar+Ag{1}*Rdif{1}+(Bg{1}+Cg{1})*Rdif{1}*Pbar+(Bg{1}+Cg{1})*Rbar*Pdif{1});
n2 =norm(Ag{2}*Sbar1+Vg{2}+Ag{1}*S1dif{1});
n3 =norm(Ag{2}*Sbar2+Ag{1}*S2dif{1});
n4 =norm(Bg{2}*(Rbar*Qbar1+Sbar1)+Wg{2}+Bg{1}*Rdif{1}*Qbar1+Bg{1}*(Rbar*Q1dif{1}+S1dif{1}));
n5 =norm(Bg{2}*(Rbar*Qbar2+Sbar2)+Bg{1}*Rdif{1}*Qbar2+Bg{1}*(Rbar*Q2dif{1}+S2dif{1}));
n = norm([n1 n2 n3 n4 n5]);
disp(['Norm of residuals from equations describing differentials ' num2str(n) ]);