%% Model selection and Out-of-Sample $R^2$
% This exercise looks at the effect of model selection on out-of-sample $R^2$.  
% The true model has 3 variables $y_t = \beta_1x_{1,t-1}+\beta_2x_{2,t-1}+\beta_3x_{3,t-1}+\epsilon_t$.  
% These three variables are correlated with a correlation of $\rho$.  In addition, 
% there are 5 more variables that are available which are correlated with the 
% true variables with a correlation of $\rho$.  When $\rho$ is large, these variables 
% will look similar to the true variables.  The ${\beta}$ is generated with two 
% considerations in mind:
% 
% * There is a parameter $\lambda\in(0,1]$ (called strength) that determines 
% how the coefficients decline in importance.  The coefficient are constructed 
% so that $\beta_j/\beta_{j-1} = \lambda$. When $\lambda=1$ all ${\beta}$ will 
% have the same value.  
% * The variance of residual is set of 1.  The ${\beta}$ are scaled so that 
% $R^2 =Var(x{\beta}) / (Var(x{\beta}) + Var(\epsilon))$ where $R^2$ is a model 
% parameter. 
% 
% All data in this exercise is simulated. We will conduct 8 scenarios using 
% all permutations of:
% 
% * Correlation $\rho\in\{0.2,0.8\}$
% * Strength $\lambda\in\{0.6,1\}$
% * $R^2\in\{0.1,0.25}$
% * 
%% Model Selection
% I will consider 5 model selection methods:
% 
% * AIC
% * BIC
% * 10-fold Cross Validation
% * GtS
% * StG
% 
% 
% 
% When using AIC, BIC and X-validation, all possible model with the regressors 
% will be evaluated leading to a total of $2^8-1$ models since I won't consider 
% models with no regressors. 
% 
% First I restore the random number generator to its default so that multiple 
% runs of the simulation will generate similar answers.

rng('default')
%% 
% Next I setup the selection arrays which will contains the entire set of 
% possible column combinations that could be jointly setup.


selection = {};
for i = 1:8
    vals = nchoosek(1:8,i);
    vals = mat2cell(vals,ones(size(vals,1),1),i);
    selection = [selection;vals];
end
%% Single Run
% Before writing loops to compute all simulations, it is useful to step through 
% the steps for a single run. The data is generated using the function |generate_data|. 
% The model is selected and fit using half and then will be evaluated using the 
% second half.
% 
% *Note*: See the functions for information on the steps used to implement 
% the various statistical measures
%%
n = 250;
r2 = 0.25;
rho = 0.8;
strength=0.6;
[y, x] = generate_data(n, 3, 5, r2, strength, rho);
full_y = y;
full_x = x;
y = y(1:n/2);
x = x(1:n/2,:);
%% Information Criteria
% The information criteria are computed in a helper function.  Here I compute 
% the IC and show the columns selected.  Ideally these would both be [1, 2, 3] 
% since these are the true regressors.

[aic, bic] = compute_ic(y,x,selection);
aic_cols = selection{aic == min(aic)}
bic_cols = selection{bic == min(bic)}
%% Cross validation
% Cross varidation is also implemented in a function.

folds = 10;
sse = kfold_cross_val(y,x,folds,selection);
xval_cols = selection{sse == min(sse)}

%% General-to-Specific and Specific-to-General
% Finally GtS and StG are implemented in their own functions. 

gts_cols = gts(y,x,0.05)
stg_cols = stg(y, x, 0.05)
%% Out-of-Sample Measure
% The final step is to use the columns selected by the different methods to 
% compute in-sample $\beta$ and then to evaluate these in-sample $\beta$ using 
% the out-of-sample data. 
% 
% I use a table to format output before showing it. 
%%
cols = {aic_cols, bic_cols, xval_cols, gts_cols, stg_cols};
oos_sse = zeros(5,1);
is_sse = zeros(5,1);
% Split the sample into eval and fit
eval_y = full_y(n/2+1:n) ;
eval_x = full_x(n/2+1:n, :);
fit_y = y(1:n/2);
fit_x = x(1:n/2, :);
% Loop across methods
for i = 1:5
    b = x(1:n/2, cols{i}) \ y(1:n/2);
    is_e = fit_y - fit_x(:, cols{i}) * b;
    is_sse(i) = is_e'*is_e;
    e = eval_y - eval_x(:, cols{i}) * b;
    oos_sse(i) = e'*e;
end
oos_r2 = 1 - oos_sse ./ (eval_y' * eval_y);
is_r2 = 1 - is_sse ./ (fit_y' * fit_y);
% Use a table to display results
names = {'AIC','BIC','X-Validation','GtS','StG'};
t = table(oos_sse, oos_r2, is_r2,  'RowNames',names,...
        'VariableNames',{'OutOfSample_SSE', 'OutOfSample_R2', 'InSample_R2'});
disp(t)
%% Full run
% We can now wrap the same code in a loop to execute across multiple configs. 
% 
% First I setup the configs for different 
%%
configs = {};
for strength = [0.5, 1]
    for rho = [0.2, 0.8]
        for r2 = [0.1, 0.2]
            configs = [configs;[strength, rho, r2]];
        end
    end
end
%% 
% The final step is to loop over the configs to run the models. |B| is the 
% number of runs of the simulation.  I keep it small so that this will run quickly 
% enough. 
% 
% The results are printed using a table.  I report three measures:
% 
% * The average in-sample SSE 
% * The average out-of-sample SSE
% * The average in-sample R2
% * The average out-of-sample R2
% * The average number of variables selected
% * The percentage of models that selected only the three relevant variables

B = 500;
for ii = 1:length(configs)
    c = configs{ii};
    strength = c(1);
    rho = c(2);
    r2 = c(3);
    
    oos_r2 = zeros(5,B);
    is_r2 = zeros(5,B);
    oos_sse = zeros(5,B);
    is_sse = zeros(5,B);
    num_var = zeros(5,B);
    correct_sel = zeros(5,B);
    for bb = 1:B
        [y, x] = generate_data(n, 3, 5, r2, strength, rho);
        full_y = y;
        full_x = x;
        fit_y = y(1:n/2);
        fit_x = x(1:n/2,:);
        eval_y = y(n/2+1:n);
        eval_x = x(n/2+1:n,:);        
        [aic, bic] = compute_ic(y,x,selection);
        aic_cols = selection{aic == min(aic)};
        bic_cols = selection{bic == min(bic)};
        gts_cols = gts(y,x,0.05);
        
        sse = kfold_cross_val(y,x,folds,selection);
        xval_cols = selection{sse == min(sse)};
        stg_cols = stg(y, x, 0.05);
        cols = {aic_cols, bic_cols, xval_cols, gts_cols, stg_cols};
        
        for i = 1:5
            b = fit_x(:, cols{i}) \ fit_y;
            e = fit_y - fit_x(:, cols{i}) * b;
            is_sse(i,bb) = e'*e;
            is_r2(i,bb) = 1 - is_sse(i,bb) / (fit_y'*fit_y);
            
            e = eval_y - eval_x(:, cols{i}) * b;
            oos_sse(i,bb) = e'*e;
            oos_r2(i,bb) = 1 - oos_sse(i,bb) / (eval_y'*eval_y);
            num_var(i,bb) = length(cols{i});
            correct_sel(i,bb) = isempty(setxor(cols{i},1:3));
        end
    end
    mean_sse = mean(oos_sse,2);
    disp('-------------------------------------------------------------')
    t = table(c(1),c(2), c(3), 'VariableNames',{'Strength', 'Rho', 'R2'});
    disp(t)
    disp('')
    disp(['Best selection is: ' names{mean_sse == min(mean_sse)}])
    disp('')
    t = table(mean_sse, mean(oos_r2,2), mean(is_sse,2), mean(is_r2,2), mean(num_var,2),mean(correct_sel,2),...
        'RowNames',names,...
        'VariableNames',{'OutOfSample_SSE', 'OutOfSample_R2', 'InSample_SSE', 'InSample_R2', 'NumVar', 'PrCorrect'});
    t = sortrows(t);
    disp(t)
end
%% Conclusions
% Overall these methods work well and the out-of-sample SSE is close to the 
% true one (125).  When 
% 
% Putting the results together, there are a few patterns that emerge:
% 
% * In-sample $R^2$ is always larger than out-of-sample $R^2$ , often substantially. 
% * AIC and cross validation tend to select larger models
% * When the regressor correlation is low then it is easier to get the correct 
% model.  
% * When strength is high (all $\beta$  are equal) ...
% * Overall the better selectors in this exercise were the looser ones, cross 
% validation and GtS.
% * When correlation is high, small models can work well since some of the information 
% in other regressors are in the small number of included regressors.
% 
% *Note* these is the dependent on this setup.  Other setups will provide 
% other rankings. This is in particular a simple exercise where the data are all 
% IID.  Adding dependence and heteroskedasticity may change these rankings.