%{
    run_all.m
    ---------

    This code estimates the income process for states using the Kalman
    Filter.

    The model is a random walk permanent plus iid transitory process.

    Model: 
            y_t  = rw_t + u_t
            rw_t = rw_{t-1} + e_t 
              
    State-Space Representation:
            x_t = K x_{t-1} + R v_t

            v_t = (e_t; u_t)
            x_t = (e_t; u_t-u_{t-1}; - u_{t})

            y_t = Z x_t 

                (0 0 0)                 (1  0)              (1)
            K = (0 0 1)             R = (0  1)          Z = (1)
            	(0 0 0)                 (0 -1)              (0)

            Q = (sigma1    0   )
                (   0    sigma2)

%}

clc
clear
close all

%% Set up project directory.
%{
    if ismac
        id = getenv('USER');
        main = strcat('/Users/',id,'/Dropbox/Intra_Macro/');                % Project Folder. Change to your folder.
        cd (strcat(main))
    else
        id = getenv('USERNAME');
        main = strcat('/Users/',id,'/Dropbox/Intra_Macro/');                % Project Folder. Change to your folder.
        cd (strcat(main))
    end

    addpath(genpath(main))                                                  % Add all subfolders to path. Note: Make sure file names are unique.
    estout = strcat(main,'Output/Estimates/Regression/');                   % Folder for diary file.
    datain = strcat(main,'Data/Final/Excel/ydata.xlsx');                    % Variables file from Stata.
    parsout = strcat(main,'Data/Final/Excel/pars.xlsx');                    % Export estimates.
    pincout = strcat(main,'Data/Final/Excel/pdata.xlsx');                   % Export permanent income.
%}    
%% Load and clean data.

    dat = kalman.cleaner(datain);                                           % Replace 'datain' with the data file (and the path if in different folder).

%% Use the Kalman filter.

    dat.pinit = 0;                                                          % Initial estimation.
    x0 = [0.5 0.5];
    
    options = optimset('Display','iter');
    kfilter0 = @(x)(kalman.filter(x,dat));
    xhat1 = fminunc(kfilter0,x0,options);                                   % Estimates.
    
    dat.pinit = 1;                                                          % Use stationary distribution.
    kfilter = @(x)(kalman.filter(x,dat));

    options = optimset('Display','iter');
    [sigmahat,~,~,~,~,H] = fminunc(kfilter,abs(xhat1),options);             % Estimates.
    sigmahat = abs(sigmahat);
    
    clear x0 xhat1

%% Use the Kalman Smoother.
   
    clc
    dat = kalman.smoother(sigmahat,dat);
    
%% Standard errors.

    % Hessian S.E. 
    
    se_hess = diag(inv(H)).^(1/2);                                          % S.E. from Empirical Hessian.
            
    dat = rmfield(dat,'pinit');   
    
%% Output.
   
	clc
    
	%diary([estout,'Kalman_',datestr(now,'dd-mm-yy','local'),'_',datestr(now,'hh-MM-ss','local'),'.txt'])
    
    diary HW1
    
    fprintf('----------------------------------------------------------\n')
    fprintf('Kalman Filter Estimates:\n\n')
    
    fprintf('sigmap \t %0.4f \n',sigmahat(1))
    fprintf('SE \t\t %0.4f \n\n',se_hess(1))
    
    fprintf('sigmat \t %0.4f \n',sigmahat(2))
    fprintf('SE \t\t %0.4f \n',se_hess(2))
    
	diary off
 
%% Extract variables from matrix.

    dat.pshocks = zeros(dat.T,dat.N);
    
    for i = 1:dat.N                                                         % Permanent shocks.
        dat.pshocks(:,i) = dat.alsmsave(:,3*(i-1)+1);
    end

    dat.tshocks = zeros(dat.T,dat.N);
    
    for i = 1:dat.N                                                         % Delta transitory shocks.
        dat.dtshocks(:,i) = dat.alsmsave(:,3*(i-1)+2);
    end

    dat.tshocks = zeros(dat.T,dat.N);
    
    for i = 1:dat.N                                                         % Transitory shocks.
        dat.tshocks(:,i) = dat.alsmsave(:,3*(i-1)+3);
    end
    
    clear dat.alsmsave
   
%% Generate permanent income.

    perm = zeros(dat.T,dat.N);
    p0   = dat.y0;
    perm(1,:) = dat.g + dat.pshocks(1,:) + p0;
    
    for i = 2:dat.T
        perm(i,:) = dat.g + dat.pshocks(i,:) + perm(i-1,:);
    end
    
    dat.perm = perm;    
    dat.perm = exp(reshape(dat.perm,[],1));
    
    clear perm
       
%% Export data. 

    dat.panelid = reshape(dat.panelid,dat.T+1,dat.N);
    dat.panelid = dat.panelid(2:end,:);
    dat.panelid = reshape(dat.panelid,[],1);
    
    dat.timeid = reshape(dat.timeid,dat.T+1,dat.N);
    dat.timeid = dat.timeid(2:end,:);
    dat.timeid = reshape(dat.timeid,[],1);
            
    A = {'panelid', 'timeid', 'perm'};
    B = num2cell(dat.panelid,3);
    C = num2cell(dat.timeid,3);
    D = num2cell(dat.perm,3);
    
    if exist(pincout, 'file')
        delete(pincout);
    end
    
    writecell([A;B C D],pincout)
    
%% Export estimates.    
    
    A = {'sigmap', 'sigmat', 'G'};
    B = num2cell(abs(sigmahat),1);
    C = mean(exp(dat.g));
    
    if exist(parsout, 'file')
        delete(parsout);
    end

    writecell([A;[B C]],parsout)
    
    clear A B C D