API tools faq
paste
Advertisement
thanhqtran

RBC basic

thanhqtran
Jul 3rd, 2025 (edited)
136
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
MatLab 4.39 KB | Software | 0 0
raw download clone embed print report
  1. %rbc uhlig
  2. clear all;
  3.  
  4. %params
  5. alpha = 0.35;
  6. beta = 0.985;
  7. eta = 2;
  8. phi = 1.5;
  9. delta = 0.025;
  10. rhoa = 0.95;
  11.  
  12. % Step 1: Compute the steady-state of all variables
  13. rss = 1 / beta + delta - 1;
  14. wss = (1-alpha) * ((alpha/rss)^(alpha/(1-alpha)));
  15. nss = ((wss^(1/eta))/( wss/(1-alpha) - delta*((wss/(1-alpha))^(1/alpha) ) ))^(1/(phi/eta + 1));
  16. iss = delta*((alpha/rss)^(1/(1-alpha)))*nss;
  17. css = (wss/(nss^phi))^(1/eta);
  18. kss = ((alpha/rss)^(1/(1-alpha)))*nss;
  19. yss = kss^alpha * nss^(1-alpha);
  20. disp([rss, wss, yss, iss, css, nss, kss]);
  21.  
  22.  
  23. % Define matrices
  24. A = [1; 0; 0; 0; 0; 0];
  25. B = [-(1-delta); -alpha; 0; 0; 1; 0];
  26. C = [0, 0, 0, 0, 0, -delta; ...
  27.      1, 0, -(1-alpha), 0, 0, 0; ...
  28.      0, eta, phi, -1, 0, 0; ...
  29.      -1, 0, 1, 1, 0, 0; ...
  30.      -1, 0, 0, 0, 1, 0; ...
  31.      -yss, css, 0, 0, 0, iss];
  32. D = [0; -1; 0; 0; 0; 0];
  33. F = [0]; G = [0]; H = [0];
  34. J = [0, eta/beta, 0, 0, -rss, 0];
  35. K = [0, -eta/beta, 0, 0, 0, 0];
  36. L = [0];
  37. M = [0];
  38. N = [rhoa];
  39.  
  40. % Solve for P (quadratic equation)
  41. C_inv = inv(C);
  42. a = F - J * C_inv*A;
  43. b = - (J * C_inv * B - G + K * C_inv * A);
  44. c = - K * C_inv * B + H;
  45. DELTA = b^2 - 4 * a * c;
  46. P1 = (-b + sqrt(DELTA))/(2*a);
  47. P2 = (-b - sqrt(DELTA))/(2*a);
  48. P = min(abs(P1), abs(P2));
  49.  
  50. % Solve for R
  51. R = -C_inv * (A * P + B);
  52.  
  53. % Solve for Q
  54.  
  55. k = 1; % Number of columns in Q, based on the dimension of z_t
  56. I_k = eye(k);
  57. % LHS matrix for the system
  58. LHS = kron(N', F - J * C_inv * A) + kron(I_k, J * R + F * P + G - K * C_inv * A);
  59. % RHS vector for the system
  60. RHS = (J * C_inv * D - L) * N + K * C_inv * D - M;
  61. % Solve for vectorized Q
  62. Q_vec = inv(LHS) * vec(RHS);
  63. %Q_vec = RHS/LHS;
  64. % Since Q is 1x1, assign directly
  65. Q = Q_vec; % No need to reshape for scalar values
  66.  
  67. % Solve for S
  68. S = -C_inv * (A * Q + D);
  69.  
  70. % Display results
  71. disp('P = '); disp(P);
  72. disp('R = '); disp(R);
  73. disp('Q = '); disp(Q);
  74. disp('S = '); disp(S);
  75.  
  76. % Run uhlig_solver.m first
  77. % Impulse Response
  78. rhoa = .95;
  79.  
  80. % one-time shock with epsilon = .01
  81. T = 80; %time horizon
  82.  
  83. tilde_k = zeros(1,T);
  84. tilde_y = zeros(1,T);
  85. tilde_c = zeros(1,T);
  86. tilde_n = zeros(1,T);
  87. tilde_w = zeros(1,T);
  88. tilde_r = zeros(1,T);
  89. tilde_i = zeros(1,T);
  90. tilde_A = zeros(1,T);
  91.  
  92. %% IRF
  93. % one-time shock
  94. tilde_A(1) = 0.01; %standard error of shock
  95.  
  96. for t = 1:T-1
  97.     tilde_A(t+1) = rhoa * tilde_A(t);
  98.     tilde_k(t+1) = P * tilde_k(t) + Q * tilde_A(t);
  99.     tilde_y(t) = R(1) * tilde_k(t) + S(1) * tilde_A(t);
  100.     tilde_c(t) = R(2) * tilde_k(t) + S(2) * tilde_A(t);
  101.     tilde_n(t) = R(3) * tilde_k(t) + S(3) * tilde_A(t);
  102.     tilde_w(t) = R(4) * tilde_k(t) + S(4) * tilde_A(t);
  103.     tilde_r(t) = R(5) * tilde_k(t) + S(5) * tilde_A(t);
  104.     tilde_i(t) = R(6) * tilde_k(t) + S(6) * tilde_A(t);
  105. end
  106.  
  107. variables = {tilde_y,tilde_c,tilde_n,tilde_w,tilde_r,tilde_i, tilde_k, tilde_A};
  108. labels = {'y', 'c', 'n', 'w', 'r', 'i', 'k', 'A'};
  109. horizon = s+1:T-1; %choose the horizon after the shock
  110.  
  111. figure;
  112. for i = 1:length(variables)
  113.     subplot(3,3,i);
  114.     plot(horizon, variables{i}(horizon), 'b', 'LineWidth', 1); %
  115.     hold on;
  116.     title(labels{i});
  117.     yline(0, 'Color','r', 'LineWidth', 1); % Zero line in red
  118.     grid on;
  119. end
  120.  
  121. %% Stochastic simulation
  122. rhoa = .95;
  123. sigmae = 0.01;
  124.  
  125. % one-time shock with epsilon = .01
  126. T = 1000; %nbumber of periods
  127.  
  128. tilde_k = zeros(1,T);
  129. tilde_y = zeros(1,T);
  130. tilde_c = zeros(1,T);
  131. tilde_n = zeros(1,T);
  132. tilde_w = zeros(1,T);
  133. tilde_r = zeros(1,T);
  134. tilde_i = zeros(1,T);
  135. tilde_A = zeros(1,T);
  136.  
  137. % a series of random shocks
  138. randn('seed',666);
  139. e = sigmae*randn(1,T);
  140. % the time series
  141. for t = 1:T-1
  142.     tilde_A(t+1) = rhoa * tilde_A(t) + e(t);
  143.     tilde_k(t+1) = P * tilde_k(t) + Q * tilde_A(t);
  144.     tilde_y(t) = R(1) * tilde_k(t) + S(1) * tilde_A(t);
  145.     tilde_c(t) = R(2) * tilde_k(t) + S(2) * tilde_A(t);
  146.     tilde_n(t) = R(3) * tilde_k(t) + S(3) * tilde_A(t);
  147.     tilde_w(t) = R(4) * tilde_k(t) + S(4) * tilde_A(t);
  148.     tilde_r(t) = R(5) * tilde_k(t) + S(5) * tilde_A(t);
  149.     tilde_i(t) = R(6) * tilde_k(t) + S(6) * tilde_A(t);
  150. end
  151.  
  152. variables = {tilde_y,tilde_c,tilde_n,tilde_w,tilde_r,tilde_i, tilde_k, tilde_A, e};
  153. labels = {'y', 'c', 'n', 'w', 'r', 'i', 'k', 'A', 'e'};
  154. horizon = s:200; %choose the horizon after the shock
  155.  
  156. figure;
  157. for i = 1:length(variables)
  158.     subplot(3,3,i);
  159.     plot(horizon, variables{i}(horizon), 'b', 'LineWidth', 1); %
  160.     hold on;
  161.     title(labels{i});
  162.     yline(0, 'Color','r', 'LineWidth', 1); % Zero line in red
  163.     grid on;
  164. end
Tags: uhlig
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!