% Monte Calro Simulations for Choi and Jeong (2019)
% Last checking: 05/20/2019

clear all
rng('default');% random number generator

% Data environment       
kmax=8;% Maximum number of factors
T=150;% # of time series observations
N=150;% # of crosssectional observations
r=5;% # of true factors
h=12;% forecasing horizon
nc=0.9;N1=floor(N*nc);% # of I(1) observations, nc=N1/N

% True parameters
alpha=1;% Autoregressive coefficient of Factor, set to be 1, i.e., F_t=F_t-1+u_t
sigmalam=1;% Variance of the True Factor Loading (for \lambda_i)
sigmalam_s=1; % Variance of the true factor loading for stationary part (for \lambda_i^s)
Lam=sqrt(sigmalam)*randn(N1,r);% Generating the true factor loading matrix
Lam_s=sqrt(sigmalam_s)*randn(N-N1,r);
iter=2000;% # of iterations (#=2,000)
ferr_Level=zeros(iter,1);ferr_Differenced=zeros(iter,1);% bins

for jj=1:iter

% True Model (r,sigmalam,alpha,beta,t,n,y0)
throw=60;% We generate T+throw(=60) observations for each i and employed the last T observations for each series
[F, F1, F0]=genar(alpha,T,r,throw);% Generating the true factor matrix;
CD=F*(Lam');% True common component component (Nonstationary part)
CD_s=(F-F1)*(Lam_s');% True common component (Stationary part)

% Error term specification
e=generroreffi(T,N1,r,throw);% idiosyncratic part for the nonstationary part
%e=genar(1,T,N1,throw);% unit root idiosyncratic errors
e_s=generroreffi(T,N-N1,r,throw);% idiosyncratic part for the stationary part
X=CD+e;% True factor model (Nonstationary part)
x_s=CD_s+e_s;% True factor model (Stationary part)
DX=[X(2:T,:)-X(1:T-1,:) x_s(2:T,:)];% Differenced data (T-1) by N (N=N1+N0)
T_1=T-1;% # of time series obervations for the differenced data
II=ones(T,1);% for the constant regressor

% Generating the true forecasting equation: y_t+h=\mu+a*F_t+b*y_t+epsil_t+h
y1=genar(1,h,1,throw);% initial h values for the forecasting equation
b=1;a=ones(r,1);mu=1;% parameters in the forecasting equation
epsil=randn(T,1);% errors in the forecasting equation
%epsil=genar(0.5,T,1,throw);% AR(1) errors in the forecasting equation
y=zeros(T,1);% bin for y_t+h
for tt=1:h
y(tt)=y1(tt); % for the first h observations, y_1, ... y_h, they are generated by the unit root process
end
for tt=1+h:T
y(tt)=mu+y(tt-h)*b+F(tt-h,:)*a+epsil(tt);
end
yh=y(h+1:end);% level y
dyh=yh(2:end)-yh(1:end-1);% differenced y

% Estimation of factor spaces
  fhat=estimation_fac(X(1:end-h,:),r,0);% estimated I(1) factor spaces
  Dfhat=estimation_fac(DX(1:end-h,:),r,1);% estimated I(0) factor spaces
  
% Forecasting regression
  RE_Level=[ones(length(y(1:end-h)),1) y(1:end-h) fhat]; % regressors (level)
  RE_Differenced=[ones(length(y(2:end-h)-y(1:end-h-1)),1) y(2:end-h)-y(1:end-h-1) Dfhat]; % regressor (differenced)
  delta_Level=(RE_Level'*RE_Level)\RE_Level'*yh;% regression coefficients (level)
  delta_Differenced=(RE_Differenced'*RE_Differenced)\RE_Differenced'*dyh;% regression coefficients (differenced)
  Fitted_Level=RE_Level*delta_Level;Fitted_Differenced=RE_Differenced*delta_Differenced;% fitted values
  ferr_Level(jj)=yh(end)-Fitted_Level(end);ferr_Differenced(jj)=yh(end)-(yh(end-h)+sum(Fitted_Differenced(end-h+1:end)));% forecasting errors

end

% Results
RMSE_Level=sqrt((ferr_Level'*ferr_Level)/iter);% root mean square forecasting errors (level)
RMSE_Differenced=sqrt((ferr_Differenced'*ferr_Differenced)/iter);% root mean square forecasting errors (differenced)
OUTPUT=[RMSE_Level, RMSE_Level/RMSE_Differenced];

save('simulation.mat')




