value of type 'sym' is not convertible to 'double' (2024)

%% COMBINED MODEL

clear; close all; clc;

%% Parameters

rho = 0.04;

theta = 1.5;

nu = 0.35;

delta = 0.05;

A=1;

L=1;

alpha = 0.5;

gamma = 0.5;

n = 0.01;

%% Functions

syms kc

Ec = (A*kc.^gamma*L.^(1-gamma)) ;

syms kp

Ep = (A*kp.^gamma*L.^(1-gamma)) ;

E = Ec + Ep;

fs = A*(E.^gamma)*(L.^(1-gamma));

f_p1s_c = diff(fs, kc);

f_p1s_p = diff(fs, kp);

f = matlabFunction(fs,'vars',{kc kp});

f_p1_c = matlabFunction(f_p1s_c,'vars',{kc kp});

f_p1_p = matlabFunction(f_p1s_p,'vars',{kc kp});

u = @(c) (integral(exp(-(p-n))*c.^(1-theta)/(1-theta)));0;inf ;

% %% Steady state

% % SS for Ramsey, GR for Solow

% % Numerical solution using unidimensional Newton-Raphson

%

% metric=1; i=0; kSS=1;

% while (metric>=1e-10)

% kSS = ((delta +rho)/(A*alpha)).^(1/alpha - 1) ;

% metric = ((delta +rho)/(A*alpha)).^(1/alpha - 1); i=i+1;

% end

% cSS = f(kSS) - (delta+n)*kSS;

kSS_c = ((delta +rho)/(A*gamma)).^(1/(gamma - 1)) ;

kSS_p = ((delta +rho)/(A*gamma)).^(1/(gamma - 1)) ;

%kSS = [kSS_c kSS_p]

cSS = f(kSS_c, kSS_p) - (delta+n)*kSS_c ;

%% Grid, differentiation operators

Ik = 2000; kmin = 1e-5*kSS_c; kmax = 5*kSS_c;

k = linspace(kmin,kmax,Ik)';

dk = (kmax-kmin)/Ik;

d1Mx = spdiags(repmat([-1 sparse(1,1) 1],Ik,1),-1:1,Ik,Ik)/2;

d1Mx(1,1)=-1; d1Mx(1,2)=1; d1Mx(Ik,Ik-1)=-1; d1Mx(Ik,Ik)=1;

d1Mx = d1Mx/dk;

%% Main loop - Euler equation + Newton-Raphson (modified reflection method)

% Initialisation

metric=1; i=0;

c=k;

% Static (c-independent) component of the Jacobian

Js = spdiag(f(kc, kp) - (delta+n)*kc)*d1Mx*theta - ...

spdiag(f_p1_c(kc, kp)-delta - rho);

Jd = sparse(Ik,Ik);

tic

while metric > 1e-10

% Dynamic (c-dependent) component of the Jacobian

% Jd = spdiags([-0.5*c [0; 0.5*(c(1:N-2)-c(3:N)); 0] 0.5*c],-1:1,N,N)';

% Jd(1,1) = 2*c(1)-c(2); Jd(1,2) = -c(1);

% Jd(N,N-1) = c(N); Jd(N,N) = c(N-1) - 2*c(N);

Jd = -(spdiag(d1Mx*c)+spdiag(c)*d1Mx)*theta;

% Merit function/Error vector (from the Euler equation)

err = (d1Mx*c).*(f(k) - (delta+n)*k - c)*theta - ...

(f_p1(k)-delta - rho).*c;

c = c - (Js+Jd)\err;

metric = max(abs(err));

i=i+1;

end

toc

if metric<=1e-10

fprintf('Solver converged in %d iterations.\n',i);

else

fprintf('Solver failed to converge in %d iterations!\n',i);

return;

end

s = f(k) - (delta+n)*k - c;

%% Simple sparse diagonalisation

function[x] = spdiag(v)

[rv,cv] = size(v);

if (cv==1)

x = spdiags(v,0,rv,rv);

return;

elseif (rv==1)

x = spdiags(v',0,cv,cv);

return;

else x=v; return;

end

end

Ran in: