dataset = readtable('Master Thesis dataset (1).xlsx', 'Sheet', 'Dataset','PreserveVariableNames', true);
data = table2array(dataset(:,4:20));

T_upr = length(data);
n = size(data,2);

% Data indexing

isi         = log(1+data(:,1));        %Investor sentiment index
age         = log(1+data(:,2));      % Company age
UPR         = data(:,3);             % Underpricing rate
employees   = data(:,4);         % # of Employees
issue_size  = data(:,5);        % Issue/offer size in millions
issue_price = log(1+data(:,6));   % Issue/offer price of IPO    
MAUP        = data(:,7);        %Market-adjusted underpricing rate
Dummy       = data(:,9);        % Hot/Cold industry dummy

% Six industry dummies; Energy, industrial, consumer cyclical and
% non-cyclical, technology:
ind1        = data(:,10);
ind2        = data(:,11);
ind3        = data(:,12);
ind4        = data(:,13);
ind5        = data(:,14);
ind6        = data(:,15);

% Data/variables for regression with upr as dependent variable and isi
% as only independent variable
data_all_isi = [isi, MAUP];
data_all_isi = array2table(data_all_isi);
data_all_isi.Properties.VariableNames = {'isi','MAUP'};

% Data/variables for regression with upr as dependent variable 
data_all_upr = [isi, age, employees, issue_size, issue_price, Dummy, ind1,ind2,...
   ind3,ind4,ind5, ind6, UPR];
data_all_upr = array2table(data_all_upr);
data_all_upr.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size', 'issue_price', 'Dummy',...
   'ind1','ind2','ind3','ind4','ind5','ind6' 'upr'};

% Data/variables for regression with maup as dependent variable 
data_all_maup = [isi, age, employees, issue_size, issue_price, Dummy, ind1,ind2,...
   ind3,ind4,ind5, ind6, MAUP];
data_all_maup = array2table(data_all_maup);
data_all_maup.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size', 'issue_price', 'Dummy',...
   'ind1','ind2','ind3','ind4','ind5', 'ind6', 'maup'};


% Regression Models

% Regression: maup as dependent variable, isi as only ind.variable
regression_isi = fitlm(data_all_isi) 
reg_u_isi = regression_upr.Residuals.Raw;
reg_u2_isi = reg_u_isi.^2;

%Regression: upr as dependent variable and all ind.variables included
regression_upr = fitlm(data_all_upr) 
reg_u_upr = regression_upr.Residuals.Raw;
reg_u2_upr = reg_u_upr.^2;

%Regression: maup as dependent variable and all ind.variables included
regression_maup = fitlm(data_all_maup) 
reg_u_maup = regression_maup.Residuals.Raw;
reg_u2_maup = reg_u_maup.^2;

%Winsorization of MAUP at 1 and 99 percentile
maup_winsor = winsor(MAUP,[1 99]);
data_all_winsor = [isi, age, employees, issue_size, issue_price, Dummy, ind1,ind2,...
   ind3,ind4,ind5, ind6, maup_winsor];
data_all_winsor = array2table(data_all_winsor);
data_all_winsor.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size', 'issue_price', 'Dummy',...
   'ind1','ind2','ind3','ind4','ind5', 'ind6', 'maup'};

regression_winsor = fitlm(data_all_winsor) 
reg_u_winsor = regression_winsor.Residuals.Raw;
reg_u2_winsor = reg_u_winsor.^2;


%
% Testing for Normality within Variables

data_Normality = [isi, age, employees, issue_size, issue_price, MAUP];


N_MAUP          = jbtest(MAUP)
N_isi           = jbtest(isi)
N_age           = jbtest(age)
N_employees     = jbtest(employees)
N_size          = jbtest(issue_size)
N_price         = jbtest(issue_price)

Norm_data = [N_MAUP, N_isi, N_age, N_employees, N_size, N_price];
tbl_norm = array2table(Norm_data)
tbl_norm.Properties.VariableNames = {'MAUP', 'isi', 'age', 'employees', 'issue_size', 'issue_price'};
tbl_norm.Properties.RowNames = {'Normality = 1, 0 otherwise'};
disp(tbl_norm)

%
% Heteroscedasticity - White's test

% Create squared level variables 
isi_2         = isi.^2;
age_2         = age.^2;
employees_2   = employees.^2;
issue_size_2  = issue_size.^2;
issue_price_2 = issue_price.^2;
upr_2         = UPR.^2;
maup_2        = MAUP.^2;


% Specify variables to be used in aux regression for White's test (upr)
aux_reg_upr_data = [isi, age, employees, issue_size, issue_price,Dummy,...
    ind1,ind2,ind3,ind4,ind5,ind6,ind7,ind8,isi_2, age_2, employees_2,...
    issue_size_2, issue_price_2, reg_u2_upr];

aux_reg_upr_data = array2table(aux_reg_upr_data);
aux_reg_upr_data.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size',...
    'issue_price', 'Dummy','ind1','ind2','ind3','ind4','ind5','ind6','ind7','ind8',...
    'isi^2', 'age^2', 'employees^2', 'issue_size^2', 'issue_price^2','u^2'};

% Run auxiliary regression used in White's test (upr)
aux_reg_upr_White = fitlm(aux_reg_upr_data)

URSS_upr     = aux_reg_upr_White.SSE;
T_upr        = aux_reg_upr_White.NumObservations;
k_upr        = aux_reg_upr_White.NumEstimatedCoefficients;
R2_upr_White = aux_reg_upr_White.Rsquared.Ordinary;
alpha        = [0.1, 0.05, 0.01];
m            = 10;

%
%

% Specify variables to be used in aux regression for White's test (maup)
aux_reg_maup_data = [isi, age, employees, issue_size, issue_price,...
    isi_2, age_2, employees_2, issue_size_2, issue_price_2, reg_u2_maup];

aux_reg_maup_data = array2table(aux_reg_maup_data);
aux_reg_maup_data.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size',...
    'issue_price','isi2', 'age2', 'employees2', 'issue_size2', 'issue_price2','u2'};

% Run auxiliary regression used in White's test (upr)
aux_reg_maup_White = fitlm(aux_reg_maup_data)

URSS_maup      = aux_reg_maup_White.SSE;
T_maup         = aux_reg_maup_White.NumObservations;
k_maup         = aux_reg_maup_White.NumEstimatedCoefficients;
R2_maup_White  = aux_reg_maup_White.Rsquared.Ordinary;

% White's test: Chi^2 test

% For upr regression
tstat_White_upr   = R2_upr_White * T_upr  %Test statistic for test
cval_White_upr    = chi2inv(1-alpha,m+8) %Critical calue
%p_White_upr       = 1 - cdf('Chisquare', tstat_White_upr, m) %p-value
reject_White_upr  = tstat_White_upr > cval_White_upr %Rejection
% H0:data is homoscedastic
% H1: data is heteroscedastic

% If reject results =
% 1: we reject H0 --> Significant heteroscedasticity
% 0: Fail to reject H0 --> data is homoscedastic

%

% For maup regression
tstat_White_maup   = R2_maup_White * T_maup  %Test statistic for test
cval_White_maup    = chi2inv(1-alpha,m) %Critical calue
%p_White_maup       = 1 - cdf('Chisquare', tstat_White_maup, m) %p-value
reject_White_maup  = tstat_White_maup > cval_White_maup %Rejection

tbl_reject_White = array2table(reject_White_maup);
tbl_reject_White.Properties.VariableNames = {'1%', '5%', '10%'};
tbl_reject_White.Properties.RowNames = {'Homoscedasticity = 0'};
disp(tbl_reject_White)
% H0:data is homoscedastic
% H1: data is heteroscedastic

% If reject results =
% 1: we reject H0 --> Significant heteroscedasticity
% 0: Fail to reject H0 --> data is homoscedastic

% Autocorrelation - Breusch-Godfrey test with ten lags

% Lagged estimated residuals for upr regression
u1  = lagmatrix(reg_u_upr, 1);
u2  = lagmatrix(reg_u_upr, 2);
u3  = lagmatrix(reg_u_upr, 3);
u4  = lagmatrix(reg_u_upr, 4);
u5  = lagmatrix(reg_u_upr, 5);
u6  = lagmatrix(reg_u_upr, 6);
u7  = lagmatrix(reg_u_upr, 7);
u8  = lagmatrix(reg_u_upr, 8);
u9  = lagmatrix(reg_u_upr, 9);
u10 = lagmatrix(reg_u_upr, 10);

% Missing value of lagged residuals are set to zero
lagged_res_data = [u1, u2, u3, u4, u5, u6, u7, u8, u9, u10];
lagged_res_data = fillmissing(lagged_res_data, 'constant', 0);

%Variables to be used for auxilary regression for Breusch-Godfrey test
BGdata_upr = [isi, age, employees, issue_size, issue_price,...
    lagged_res_data, reg_u_upr];
    
BGdata_upr = array2table(BGdata_upr);
BGdata_upr.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size',...
    'issue_price','u1', 'u2', 'u3', 'u4', 'u5','u6','u7','u8','u9','u10','u'};

% Run auxiliary regression used in Breusch-Godfrey test
aux_reg_BG_upr = fitlm(BGdata_upr)

R2_aux_BG_upr  = aux_reg_BG_upr.Rsquared.Ordinary;

%

% Lagged estimated residuals for maup regression
u1_maup   = lagmatrix(reg_u_maup, 1);
u2_maup   = lagmatrix(reg_u_maup, 2);
u3_maup   = lagmatrix(reg_u_maup, 3);
u4_maup   = lagmatrix(reg_u_maup, 4);
u5_maup   = lagmatrix(reg_u_maup, 5);
u6_maup   = lagmatrix(reg_u_maup, 6);
u7_maup   = lagmatrix(reg_u_maup, 7);
u8_maup   = lagmatrix(reg_u_maup, 8);
u9_maup   = lagmatrix(reg_u_maup, 9);
u10_maup  = lagmatrix(reg_u_maup, 10);

% Missing value of lagged residuals are set to zero
lagged_res_data_maup = [u1_maup, u2_maup, u3_maup, u4_maup, u5_maup,...
     u6_maup, u7_maup, u8_maup, u9_maup, u10_maup];
lagged_res_data_maup = fillmissing(lagged_res_data_maup, 'constant', 0);

%Variables to be used for auxilary regression for Breusch-Godfrey test
BGdata_maup = [isi, age, employees, issue_size, issue_price,...
    lagged_res_data_maup , reg_u_maup];
    
BGdata_maup = array2table(BGdata_maup);
BGdata_maup.Properties.VariableNames = {'isi', 'age', 'employees', 'issue_size',...
    'issue_price','u1', 'u2', 'u3', 'u4', 'u5','u6','u7','u8','u9','u10','u'};

% Run auxiliary regression used in Breusch-Godfrey test
aux_reg_BG_maup = fitlm(BGdata_maup)

R2_aux_BG_maup  = aux_reg_BG_maup.Rsquared.Ordinary;

% Breusch-Godfrey test with ten lags: Chi^2 test

% BG test for upr 
tstat_BG_upr = R2_aux_BG_upr * (T_upr-m)
cval_BG_upr  = chi2inv(1-alpha,m)
p_BG_upr         = 1 - cdf('Chisquare', tstat_BG_upr, m)
reject_BG_upr    = tstat_BG_upr > cval_BG_upr

% H0: there is no serial correlation
% H1: There is autocorrelation

% 1: condition is true, i.e., reject H0
% 0: condition is not true, i.e., do not reject H0

%

% BG test for maup
tstat_BG_maup    = R2_aux_BG_maup * (T_maup-m)
cval_BG_maup     = chi2inv(1-alpha,m)
p_BG_maup        = 1 - cdf('Chisquare', tstat_BG_maup, m)
reject_BG_maup   = tstat_BG_maup > cval_BG_maup

%Table of BG test
tbl_reject_BG = array2table(reject_BG_maup);
tbl_reject_BG.Properties.VariableNames = {'1%', '5%', '10%'};
tbl_reject_BG.Properties.RowNames = {'Homoscedasticity = 0'};
disp(tbl_reject_BG)


% H0: there is no serial correlation
% H1: There is autocorrelation

% 1: condition is true, i.e., reject H0
% 0: condition is not true, i.e., do not reject H0

%Joint table of White and BG
data_joint = [reject_White_maup;reject_BG_maup];
tbl_joint = array2table(data_joint);
tbl_joint.Properties.VariableNames = {'1%', '5%', '10%'};
tbl_joint.Properties.RowNames = {'White test', 'Breusch-Godfrey test'};
disp(tbl_joint)

% Multicollinearity - Correlation matirx

% Correlation matrix of variables
corr_data = [UPR, MAUP, isi, age, employees, issue_size, issue_price];

Corr_matrix = corr(corr_data);
Corr_matrix_round = round(Corr_matrix,4);


tbl_corr = array2table(Corr_matrix_round,"VariableNames",...
    ["UPR","MAUP","ISI", "Issue size", "Issue price", "Age", "Employees"],...
    "RowNames",["UPR","MAUP","ISI", "Issue size", "Issue price", "Age", "Employees"]);
disp(tbl_corr)

%
% DESCRIPTIVE STATISTICS

data_DS = (dataset(:,2:3));
listing_date = data_DS.("Listing Date");
Year = data_DS.("Year");

% MAUP and UPR plotted:

figure
subplot(2,1,1); plot(listing_date,UPR)
title("Underpricing rate")
xlabel("Year")
ylabel("UPR")
subplot(2,1,2); plot(listing_date, MAUP)
title("Market-adjusted Underpricing Rate")
xlabel("Year")
ylabel("MAUP")

subplot(2,1,1)
xlim([datetime(2005,7,29,10,32,8)...
      datetime(2018,2,17,19,15,47)])
ylim([-3.8 3.5])

subplot(2,1,2)
xlim([datetime(2005,7,29,10,32,8)...
      datetime(2018,2,17,19,15,47)])
ylim([-3.8 3.5])


% Investor sentiment index

dataset_is = readtable('Master Thesis dataset (1).xlsx', 'Sheet', 'Sentiment index','PreserveVariableNames',true);
data_is = dataset_is(2:end,1:2);

is_index = data_is.Var1;
is_date  = data_is.Var2;

figure
plot(is_date,is_index)
title("Investor Sentiment Index, Jan 2005 - Dec 2018")
xlabel("Year")
ylabel("Sentiment Index")

% Variable data

isi_tbl          = data(:,1); %Investor sentiment index
age_tbl          = data(:,2); % Company age
UPR_tbl          = data(:,3); % Underpricing rate
employees_tbl    = data(:,4); % # of Employees
issue_size_tbl   = data(:,5); % Issue/offer size in millions
issue_price_tbl  = data(:,6); % Issue/offer price of IPO
MAUP_tbl         = data(:,7); %Market-adjusted underpricing rate

data_tbl = [UPR_tbl,MAUP_tbl,isi_tbl,issue_size_tbl,issue_price_tbl,age_tbl,employees_tbl];


Means   = round(mean(data_tbl,'omitnan')',4);
Median  = round(median(data_tbl,'omitnan')',3);
Stdev   = round(std(data_tbl,'omitnan')',3);
Min     = round(min(data_tbl)',3);
Max     = round(max(data_tbl)',3);


Tbl = table(Means,Median,Stdev,Min,Max);
Tbl.Properties.VariableNames = {'Mean','Median','St.dev','Min.','Max.'};
Tbl.Properties.RowNames = {'UPR', 'MAUP','ISI', 'Issue size (Mln)', 'Issue price (NOK)', 'Age (Yrs)', 'Employees'};


disp(Tbl)


% Positive investor sentiment variable data

dataset_pos = readtable('Master Thesis dataset (1).xlsx', 'Sheet', 'Pos. inv sent','PreserveVariableNames', true);
data_pos = table2array(dataset_pos(:,5:21));
data_pos = rmmissing(data_pos);

% Data variables for positive investor sentiment periods
isi_pos           = data_pos(:,1); %Investor sentiment index
age_pos           = data_pos(:,2); % Company age
UPR_pos           = data_pos(:,3); % Underpricing rate
employees_pos     = data_pos(:,4); % # of Employees
issue_size_pos    = data_pos(:,5); % Issue/offer size in millions
issue_price_pos   = data_pos(:,6); % Issue/offer price of IPO
MAUP_pos          = data_pos(:,7); %Market-adjusted underpricing rate
Dummy_pos = data_pos(:,9);
ind1_pos = data_pos(:,10);
ind2_pos = data_pos(:,11);
ind3_pos = data_pos(:,12);
ind4_pos = data_pos(:,13);
ind5_pos = data_pos(:,14);



% Negative investor variable data

dataset_neg = readtable('Master Thesis dataset (1).xlsx', 'Sheet', 'Neg. inv sent','PreserveVariableNames', true);
data_neg = table2array(dataset_neg(:,5:21));
data_neg = rmmissing(data_neg);

% Data variables for negative investor sentiment periods
isi_neg           = data_neg(:,1); %Investor sentiment index
age_neg           = data_neg(:,2); % Company age
UPR_neg           = data_neg(:,3); % Underpricing rate
employees_neg     = data_neg(:,4); % # of Employees
issue_size_neg    = data_neg(:,5); % Issue/offer size in millions
issue_price_neg   = data_neg(:,6); % Issue/offer price of IPO
MAUP_neg          = data_neg(:,7); %Market-adjusted underpricing rate
Dummy_neg = data_neg(:,9);
ind1_neg = data_neg(:,10);
ind2_neg = data_neg(:,11);
ind3_neg = data_neg(:,12);
ind4_neg = data_neg(:,13);
ind5_neg = data_neg(:,14);



%
% Difference in mean

data_pos_tbl = [UPR_pos,MAUP_pos,isi_pos,issue_size_pos,issue_price_pos,age_pos,employees_pos];
data_neg_tbl = [UPR_neg,MAUP_neg,isi_neg,issue_size_neg,issue_price_neg,age_neg,employees_neg];

Mean_pos    = round(mean(data_pos_tbl)',3);
Mean_neg    = round(mean(data_neg_tbl)',3);
Median_pos  = round(median(data_pos_tbl)',3);
Median_neg  = round(median(data_neg_tbl)',3);
Stdev_pos   = round(std(data_pos_tbl)',3);
Stdev_neg   = round(std(data_neg_tbl)',3);
Mean_diff   = round(Mean_pos-Mean_neg,3);

%Fill in t-test for mean differences 
n_pos = length(data_pos_tbl);
n_neg = length(data_neg_tbl);
std_pos = std(data_pos_tbl);
std_neg = std(data_neg_tbl);

std_pooled = sqrt(((n_pos-1)*std_pos)+((n_neg-1)*std_neg)/(n_pos+n_neg-2))';
t_stat = round(Mean_diff./(std_pooled*sqrt(1/n_pos+1/n_neg)),3);


alpha = [0.01,0.05,0.1];
df    = n_pos + n_neg -2;
pLower= alpha/2;
pUpper= 1 - alpha/2;

tCritical = tinv([pLower, pUpper],df);
p = (1 - normcdf(t_stat))*2;


%
% Table of mean difference statistics

%Positive investor sentiment
Tbl_pos = table(Mean_pos,Median_pos,Stdev_pos);
Tbl_pos.Properties.VariableNames = {'Mean','Median','St.dev'};
Tbl_pos.Properties.RowNames = {'UPR', 'MAUP','ISI', 'Issue size', 'Issue price', 'Age', 'Employees'};
disp(Tbl_pos)

%Negative investor sentiment
Tbl_neg = table(Mean_neg,Median_neg,Stdev_neg);
Tbl_neg.Properties.VariableNames = {'Mean','Median','St.dev'};
disp(Tbl_neg)

% T-test for equality of means
Tbl_ttest = table(Mean_diff,t_stat);
Tbl_ttest.Properties.VariableNames = {'Mean Difference','t-statistic'};
disp(Tbl_ttest)

%
% Regression positive sentiment periods

% Regression Variables
isi_regpos         = log(1+isi_pos);           
age_regpos         = log(1+age_pos);         
%UPR_regpos         = UPR_pos ;          
emp_regpos         = employees_pos;   
issue_size_regpos  = issue_size_pos;    
issue_price_regpos = log(1+issue_price_pos);   
MAUP_regpos        = MAUP_pos;

% Data for regression with maup as dependent variable 
data_all_pos = [isi_regpos, age_regpos, emp_regpos,...
    issue_size_regpos, issue_price_regpos, Dummy_pos, ind1_pos,ind2_pos,...
   ind3_pos,ind4_pos,ind5_pos,ind6_pos, MAUP_regpos];
data_all_pos = array2table(data_all_pos);
data_all_pos.Properties.VariableNames = {'isi', 'age', 'employees',...
    'issue_size', 'issue_price','Dummy','ind1','ind2','ind3','ind4','ind5','ind6','maup'};

regression_pos = fitlm(data_all_pos) %Regression: maup dependent variable
reg_u_pos = regression_pos.Residuals.Raw;
reg_u2_pos = reg_u_pos.^2;

%
% Regression negative sentiment periods

%Regression variables
isi_regneg          = log(1+isi_neg);          
age_regneg          = log(1+age_neg);           
%UPR_regneg          = UPR_neg;           
emp_regneg          = employees_neg;     
issue_size_regneg   = issue_size_neg;    
issue_price_regneg  = log(1+issue_price_neg);  
MAUP_regneg         = MAUP_neg;

% Data for regression with maup as dependent variable 
data_all_neg = [isi_regneg, age_regneg, emp_regneg,...
    issue_size_regneg, issue_price_regneg, Dummy_neg, ind1_neg,ind2_neg,...
   ind3_neg,ind4_neg,ind5_neg,ind6_neg, MAUP_regneg];
data_all_neg = array2table(data_all_neg);
data_all_neg.Properties.VariableNames = {'isi', 'age', 'employees',...
    'issue_size', 'issue_price','Dummy','ind1','ind2','ind3','ind4','ind5',...
    'ind6','maup'};

regression_neg = fitlm(data_all_neg) %Regression: maup dependent variable
reg_u_neg = regression_neg.Residuals.Raw;
reg_u2_neg = reg_u_neg.^2;

%
% Plotted regressions:

plot(regression_isi)
plot(regression_maup)
plot(regression_upr)
plot(regression_winsor)
plot(regression_pos)
plot(regression_neg)

%
% Functions:

function[y,varargout] = winsor(x,p)
% WINSOR     Winsorize a vector 
% INPUTS   : x - n*1 data vector
%            p - 2*1 vector of cut-off percentiles (left, right) 
% OUTPUTS  : y - winsorized x, n*1 vector
%            i - (optional) n*1 value-replaced-indicator vector
% NOTES    : Let p1 = prctile(x,p(1)), p2 = prctile(x,p(2)). (Note
%            that PRCTILE ignores NaN values). Then 
%            if x(i) < p1, y(i) = min(x(j) | x(j) >= p1)
%            if x(i) > p2, y(i) = max(x(j) | x(j) <= p2)
% EXAMPLE  : x = rand(10,1), y = winsor(x,[10 90])
% AUTHOR   : Dimitri Shvorob, dimitri.shvorob@vanderbilt.edu, 4/15/07

if ~isvector(x)
   error('Input argument "x" must be a vector')
end  
if nargin < 2
   error('Input argument "p" is undefined')
end 
if ~isvector(p)
   error('Input argument "p" must be a vector')
end  
if length(p) ~= 2
   error('Input argument "p" must be a 2*1 vector')
end  
if p(1) < 0 || p(1) > 100
   error('Left cut-off percentile is out of [0,100] range')
end  
if p(2) < 0 || p(2) > 100
   error('Right cut-off percentile is out of [0,100] range')
end  
if p(1) > p(2)
   error('Left cut-off percentile exceeds right cut-off percentile')
end  
p = prctile(x,p);
i1 = x < p(1); v1 = min(x(~i1));
i2 = x > p(2); v2 = max(x(~i2));
y = x;
y(i1) = v1;
y(i2) = v2;
if nargout > 1
   varargout(1) = {i1 | i2};
end   
end