%% First we load dynare and get the right file directory

clear all

addpath C:/dynare/4.4.2/matlab

cd('C:\Users\ErikinWest\Documents\courses (queens)\Econ 899 - SUMMER PAPER\Christiano Course\assignment 9 forecasting')

%% Step 1: Generate some data

dynare simulate

t = 1:length(x);

%% How does the simulated data look?

plot(t,da,'blue'); hold on % productivity growing around 3% (good)
plot(t,x,'red'); hold on % output gap is mean zero random walk (good)
plot(t,tau,'green') %preferences are mean zero reverting (good)
plot(t,rstar,'blue',t,r,'red'); hold on %r and rstar fluctuate around their steady state values
plot(t,pie,'red'); hold off %as does inflation...

%% Ok so now we're ready to save the data. Let's take the last 50 observations

data = oo_.endo_simul(:,150:length(x));

%use the names to create the datapuller script
names = M_.endo_names

save('data.mat','data')

clear all

%% Step 2: Estiamte the parameters

dynare forecast

% dynare is unable to estimate all of our variables. It does a very good
% job when we only give it six to do and tell it the rest. 

%% Step 3: Include forecasts from periods 42-51

% so now we focus on forecasting. First to speed up calculation we can
% start at: mode_file=forecast_mode which tells matlab to start the
% maximization procedure at the previous mode, so calculations are sped up

dynare forecast2

%% Extracting the forecasts

% when we specify the forecast=k, we will do a k-step ahead forecast.

% suppose our data series is of length T. Then there is no difference doing
% a forecast with forecast=4 and no mention of nobs or writing that 
% nobs=[T:T] with forecast=4. 

% Now suppose we write, nobs[T-1,T], forecast=4. The first row of the
% recursive forecast matrix will contain the forecast of T, T+1, T+2, and
% T+3, using periods T-1 data. If we're only concerned about the in-sample
% fit, then we should specify, nobs[x,T-1]. But of course, this assumes we
% have a k-step ahead forecast, because column 2 of the last row of the
% recursive matrix will be for period T+1... so really nobs[x,T-k] if we
% want all of the elements in the recursive forecast matrix to be linkable
% to some specific data!

% When we specify nobs[x:y] we run a recursive forecast starting with 
% using periods [1:x] to estimate the parameters, and then forecasting
% x+1,..,x+k
% , then we use periods [1:x+1] to reestimate parameters and forecast
% x+2,x+3,..,x+l
% and we do this up until [1:y] to forecast y+1,..,y+k

% note that a recursive forecast is identical to a
% rolling forecast except that the base year doesn't change, so it's an
% 'exapnding window'
% for example...
% 1-step-ahead forecasts for:......................Rolling window.................Recursive window
% 1999M1............................................1990M1--1998M12 ...............1990M1--1998M12
% 1999M2............................................1990M2--1999M1.................1990M1--1999M1
% 1999M3............................................1990M3--1999M2.................1990M1--1999M2

% in forecast2, we're using nobs=[39:47], which means that we will have
% estimates of periods 40-48 for the OSA, 41-49 for 2SA, 42-50 for 3SA, and
% 43-51 for 4SA. 


% to extract the simple end-of-sample forecast for 8 periods: 

pie_forecast =  oo_.RecursiveForecast.Mean.pie;
r_forecast = oo_.RecursiveForecast.Mean.r;
da_forecast = oo_.RecursiveForecast.Mean.da;
%this matrix is of verttical length y-x+1 (which is the numbers of periods
%in nobs), by k, the length of the forecast step. therefore if we want to
%extract the 1-step ahead forecast, we say:

pie_1sa =  pie_forecast(:,1);
pie_2sa = pie_forecast(:,2);
r_1sa = r_forecast(:,1);
r_2sa = r_forecast(:,2);
da_1sa = da_forecast(:,1);
da_2sa = da_forecast(:,2);

%no don't forget the last term of this vector is out of sample!

% now we can see how it compares to actual data

data_puller;
x =39 ; k=4; y =47;

% plot(1:6,pie_osa(1:y-x),'blue',1:6,pie(x:y),'red')

%% Comparing OSA forecast

tt = 1:length(pie_1sa);

subplot(2,2,1);
plot(39:47,pie_1sa,'blue',39:47,pie(x:y),'red')
legend('Forecast Data','Actual Data')
title('Actual Data from Periods 39-47')

subplot(2,2,2);
plot(40:48,pie_1sa,'blue',40:48,pie(x+1:y+1),'red')
legend('Forecast Data','Actual Data')
title('Actual Data from Periods 40-48')

subplot(2,2,3);
plot(39:47,r_1sa,'blue',39:47,r(x:y),'red')
legend('Forecast Data','Actual Data')
title('Actual Data from Periods 39-47')

subplot(2,2,4);
plot(40:48,r_1sa,'blue',40:48,r(x+1:y+1),'red')
legend('Forecast Data','Actual Data')
title('Actual Data from Periods 40-48')