function [corr,bet,sig,misc ] = lrcov(x,y,q,i_corr,i_bet,i_sig,parmdir)
% This version: 11/27/2016, mww
%
% This function computes longrun covariability statistics for (x,y)
%  Input
%     x: x data (T x 1) time series
%     y: y data (T x 1) time series
%     q: number of cosine transforms used for inference. Allowable values of q are 2,3, ... , 26.   
%     i_corr: indicator variable for conducting inference about LR correlation
%     i_bet: indicator variable for conducting inference about regression coefficient of y on x (beta)
%     i_sig: indicator variable for conducting inference about sigma_y|x (standard deivation of regression error)
%     parndir: string variable that is the director where the pre-computed "sigma matrices", alfds and other parameters are stored (these are supplied with the replication materials)
%  Output
%     corr:
%       corr.posterior_median = posterior median;
%       corr.posterior_mean = posterior mean;
%       corr.posterior_cs = posterior credible set (equal-tailed) each row correponds to a different level;
%       corr.conf_interval = confidence interval each row correponds to a different level;
%       corr.level = cover probability of credible sets and confidence intervals
%       corr.r2_posterior_median = posterior median of squared correlation;
%     bet
%       bet.posterior_median = posterior median;
%       bet.posterior_mean = posterior mean;
%       bet.posterior_cs = posterior credible set (equal-tailed) each row correponds to a different level;
%       bet.conf_interval = confidence interval each row correponds to a different level;
%       bet.level = cover probability of credible sets and confidence intervals
%       bet.r2_posterior_median = posterior median of squared correlation;
%     sig
%       sig.posterior_median = posterior median;
%       sig.posterior_mean = posterior mean;
%       sig.posterior_cs = posterior credible set (equal-tailed) each row correponds to a different level;
%       sig.conf_interval = confidence interval each row correponds to a different level;
%       sig.level = cover probability of credible sets and confidence intervals
%       sig.r2_posterior_median = posterior median of squared correlation;
%     misc
%       misc.X = cosine transforms of data;
%       misc.Xproj = projections of data onto cosine terms;
%       misc.Xproj_m = projections of data onto cosine terms + constant (Xproj + sample mean)
%
%
% Carry out analysis

% Some Parameters
str_tmp = [parmdir 'level_vec']; load(str_tmp);
lfparm.level_vec = level_vec;  % Level of test in percentage points = (1 - coverage probability)
                               % Currently implement for 67% CI, 90% CI, and 95% CI 
                               % Note ... this must be read from PARM directory, because the multi-dimensional arrays there are saved in and order indicated by the entries in level_vec                                
% Initialization                                 
corr.tmp = NaN;
bet.tmp = NaN;
sig.tmp = NaN;
misc.tmp = NaN;
lfparm.q = q;

% Check for allowable values
if (q < 2) | (q > 26);
 error('q must be an integer between q=2 and q=26');
end;
% Check for missing values in data
if sum(sum(isnan([x y]))) > 0;
    error('Missing values in x and/or y not allowed \n');
end;

 % Parameters
 T = size(x,1);
   
 % Step 1: Compute cosine transforms and misc statistics;
   xp = [x y];
   X = NaN(q,2);
   Xstd = NaN(q,2);
   Xproj = NaN(T,2);
   Xproj_m = NaN(T,2);
   for i = 1:2;
     [X_mean,X_cos] = xp_compute(xp(:,i),q);
     X_norm = X_cos/sqrt(X_cos'*X_cos); 
     psi = psi_compute(T,q);
     X_proj = psi*(inv(psi'*psi))*X_cos;
     X_proj_m = X_proj + X_mean;
     X(:,i) = X_cos;
     Xstd(:,i) = X_norm;
     Xproj(:,i) = X_proj;
     Xproj_m(:,i) = X_proj_m;
   end;
   data.x = xp;             % raw data
   data.X = X;              % Cosine transform
   data.Xstd = Xstd;        % Standardized Cosine transform
   data.T = T;              % Number of observations
   % Construct Invariants
   [Xstar U] = lower_tri_invariant(data.X);
   data.Xstar = Xstar;
   data.U = U;    
   
   % Carry out analysis for given problem
      if i_corr==1
        % ----- Correlation -----------
         %fprintf('computing correlation results \n');
         lfparm.number_th0 = 28;               % Number of theta files under Ho (number of null values tested 
         lfparm.problem_label = 'varcorr';     % Problem label 
         dens_h0 = @getdens_2xscale;           % Density to use under h0
         dens_h1 = @getdens_2xscale;           % Density to use under h1
         beta_compute = @beta_compute_varcorr; % Function for computing paramter of interest for this problem
         corr_rslt = corr_compute(data,lfparm,dens_h0,dens_h1,beta_compute,parmdir);
      end  
      if i_bet==1
        % ------------------ Beta -----------
         %fprintf('computing beta results (takes a few seconds) \n');
         lfparm.number_th0 = 1;              % Number of theta files under Ho (number of null values tested 
         lfparm.problem_label = 'varbeta';   % Problem label 
         dens_h0 = @getdens_tri_beta0;       % Density to use under h0
         dens_h1 = @getdens_tri_beta1;       % Density to use under h1
         dens_w_ystar = @getdens_tri_beta0;  % Density for computing posterior
         dens_ystar = @getdens_tri_inv;      % Density of maximal invariant
         beta_compute = @beta_compute_NaN;   % Function for computing paramter of interest for this problem
         beta_rslt = b_compute(data,lfparm,dens_h0,dens_h1,dens_ystar,dens_w_ystar,beta_compute,parmdir);      
      end
       
      if i_sig==1
        % ------------------ sigma -----------
         %fprintf('computing sigma results (takes a few seconds) \n');
         lfparm.number_th0 = 1;                % Number of theta files under Ho (number of null values tested 
         lfparm.problem_label = 'varstddev';   % Problem label 
         dens_h0 = @getdens_tri_stddev0;       % Density to use under h0
         dens_h1 = @getdens_tri_stddev1;       % Density to use under h1
         dens_w_ystar = @getdens_tri_stddev0;  % Density for computing posterior
         dens_ystar = @getdens_tri_inv;        % Density of maximal invariant
         beta_compute = @beta_compute_NaN;     % Function for computing paramter of interest for this problem
         stddev_rslt = stddev_compute(data,lfparm,dens_h0,dens_h1,dens_ystar,dens_w_ystar,beta_compute,parmdir);
      end

  % Save Results
   if i_corr==1;
    corr.posterior_median = corr_rslt.posterior_median;
    corr.posterior_mean = corr_rslt.posterior_mean;
    corr.posterior_cs = corr_rslt.posterior_cs;
    corr.conf_interval = corr_rslt.ci;
    corr.level = 1-(lfparm.level_vec/100)';
    corr.r2_posterior_median = corr_rslt.r2_posterior_median;
   end;
   if i_bet==1;
    bet.posterior_median = beta_rslt.posterior_median;
    bet.posterior_mean = beta_rslt.posterior_mean;
    bet.posterior_cs = beta_rslt.posterior_cs;
    bet.conf_interval = beta_rslt.ci;  
    bet.level = 1-(lfparm.level_vec/100)';
   end;
   if i_sig==1;
    sig.posterior_median = stddev_rslt.posterior_median;
    sig.posterior_mean = stddev_rslt.posterior_mean;
    sig.posterior_cs = stddev_rslt.posterior_cs;
    sig.conf_interval = stddev_rslt.ci;  
    sig.level = 1-(lfparm.level_vec/100)';
   end;
    misc.X = data.X;
    misc.Xproj = Xproj;
    misc.Xproj_m = Xproj_m;
    
end

