function y = genYn(fname,maturities)
% Construct monthly yields based on Svensson's method in eq (22) in GSW.
% Input: 
%        file name;
%        maturities in months;
% Ouput:
%        yields: differerent columns are forward rates of different maturities

% Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 7, "Encoding", "UTF-8");
% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = ["Date", "BETA0", "BETA1", "BETA2", "BETA3", "TAU1", "TAU2"];
opts.VariableTypes = ["datetime", "double", "double", "double", "double", "double", "double"];
% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
% Specify variable properties
opts = setvaropts(opts, "Date", "InputFormat", "yyyy-MM");
% Import the data
datacell = readtable("./data/GM_sv.csv", opts);
clear opts

% get rid of data with NA
datacell(any(ismissing(datacell),2), :) = [];
date = datenum(datacell{:,1});
Y = year(date); M = month(date); D = day(date);
% Pick the data from 1972/09 to 2022/04
id_start = find(100*Y+M == 197209,1,'last');
id_end = find(100*Y+M == 202204,1,'last');
Y = Y(id_start:id_end); M = M(id_start:id_end); D = D(id_start:id_end);
data = datacell{id_start:id_end,2:end} ;
mydata = [Y M D data];
% Find the end of month observation
[~,m,~] = unique(100*mydata(:,1)+mydata(:,2),'last');
mydata = mydata(m,:);
% % find the end of quarter observation
% mydata = mydata(mod(mydata(:,2),3)==0,:);
% y.time = mydata(:,1)*100 + mydata(:,2);
%Calculate annualized forward rate 
beta0 = mydata(:,4)'; beta1 = mydata(:,5)'; beta2 = mydata(:,6)'; beta3 = mydata(:,7)'; tau1 = mydata(:,8)'; tau2 = mydata(:,9)';
n = maturities/12;
beta0 = repmat(beta0,length(n),1); beta1 = repmat(beta1,length(n),1); beta2 = repmat(beta2,length(n),1); beta3 = repmat(beta3,length(n),1); tau1 = repmat(tau1,length(n),1); tau2 = repmat(tau2,length(n),1);
n = repmat(n,1,size(mydata,1));
yield1 = beta0 + beta1.*(1-exp(-n./tau1))./n.*tau1 + beta2.*((1-exp(-n./tau1))./n.*tau1-exp(-n./tau1)) + beta3.*((1-exp(-n./tau2))./n.*tau2-exp(-n./tau2));



y = yield1';

end