trintree_calswaption( Curve, V, Period, coin )

function [ bmk_swp_premium, premium, a ] = trintree_calswaption( Curve, V, Period, coin )

% This function produces the calibrated parameters for the HW model in the extended Vasicek specification.
% The reference basket is assumed to be the ATM swaption volatility matrix (Black76 model). 
% The Volatility surface matrix V is assumed to be 10y X 10y of Expiry X Maturity.
% The model parameter vector a contains the levels of a straight line volatility function
% ranging from 0 to 20y. The mean reversion parameters is held constant across the time domain.
% 
% input
% Curve : interest rate curve object
% V : volatility matrix
% Period : frequency of payments of the underlying swaps
% coin : initial condition for the model parameters

[ bmk_swp_premium, par_rate ] = vol2par_swaption( Curve, V, Period );
bmk_swp_premium=bmk_swp_premium(:);
K=par_rate(:);
CFlowDates=cfdates(Curve.Settle,addtodate(Curve.Settle,20,'year'),Period,Curve.Basis,0)';
d=CFlowDates(Period:Period:10*Period);
FwdDates=repmat(d,10,1);
%Maturity=time2date(ExerciseDates,reshape(repmat((1:10),10,1),100,1),Curve.Compounding,Curve.Basis);
Maturity=CFlowDates((2*Period):Period:(10*Period+Period));
for i=1:9
    Maturity=[Maturity;CFlowDates((i*Period+((2*Period):Period:(10*Period+Period)))')]; %#ok<AGROW>
end

%--------------------------------------------------------
% Generate the swap portfolio cash-flow Matrix
%--------------------------------------------------------

[Q1,Q2]=cfamounts(1,Curve.Settle,[Maturity(1);Maturity],...
    Period,Curve.Basis,0,[],[],[],[Curve.Settle;FwdDates]);
[CFlowT,T]=cfport(Q1,Q2);

% the first row is fictious to generate the cash-flow times
% between the settlement date and the maturity date
CFlowT(1,:)=[];

% convert to cash-flow indicators
CFlowT(CFlowT~=0)=1;

% number of swaps in the portfolio
N=size(CFlowT,1);

% check if the fwd start date is a cash-flow date, then 
% delete the corresponding cash-flow from the CF Matirx
% otherwise include the start time in the time schedule
% in order to generate the fixing of the fwd rate in the 
% trinomial tree.

for i=1:N
    if any(T==FwdDates(i)),
        CFlowT(i,find(T==FwdDates(i),1,'first'))=0;
    else
        t=find(T>FwdDates(i),1,'first');
        T=[T(1:t-1);FwdDates(i);T(t:end)];
        CFlowT=[CFlowT(:,1:t-1),zeros(N,1),CFlowT(:,t:end)];
    end
end

%--------------------------------------------------------
% Generate the tree matching portfolio cash-flow times
%--------------------------------------------------------

% let the Tree end at the max coupon date
T_plus=addtodate(T(end),12/Period,'month');
TimeSpec=hwtimespec(Curve.Settle,[T(2:end);T_plus],Curve.Compounding);
TimeSpec.Basis=Curve.Basis;
RateSpec=Curve.toRateSpec([T(2:end);T_plus]);

% alpha can be generate outside the pricing function
[T1,T2]=meshgrid(FwdDates,T);
alpha=[zeros(1,N);diff(max(yearfrac(T1,T2,Curve.Basis),0))]'.*CFlowT;
        
%--------------------------------------------------------
% Run the portfolio calibration
%--------------------------------------------------------

a = fmincon(@swaption_premium,coin,[],[],[],[],[1e-6;1e-6;1e-6],[2;2;2],[],optimset('Display','iter','algorithm','active-set'));
    
    
    function fun = swaption_premium(a)

        VolSpec = hwvolspec(Curve.Settle, [Curve.Settle;T_plus],abs(a(1:2)),T_plus, a(3));
        Tree=hwtree(VolSpec,RateSpec,TimeSpec,'method','HW2000');

        % convert the fwd-factors tree into the fwd-rates tree
        FTree=cvtree(Tree);

        %--------------------------------------------------------
        % Begin the pricing iterations
        %--------------------------------------------------------

        [~,~,Status]=trintreeshape(Tree);
        tk=length(Status);
        SwapTree = mktrintree(tk,N,Status(1:tk),0);

            for j=tk:-1:2

                Probs=zeros([Status(j-1) Status(j)]);
                iMat=sub2ind([Status(j-1) Status(j)],repmat((1:Status(j-1))',1,3),...
                    [Tree.Connect{j-1}'-1, Tree.Connect{j-1}', Tree.Connect{j-1}'+1]);
                Probs(iMat)=Tree.Probs{j-1}';
                Bond=repmat(1./Tree.FwdTree{j-1},N,1);
                FloatLeg=(alpha(:,j-1).*CFlowT(:,j)*FTree.RateTree{j-1})';
                FixedLeg=repmat((alpha(:,j-1).*CFlowT(:,j).*K)',Status(j-1),1);

                % generate the state-matrix via backward induction
                % the first term is the expected present value of the state vector at time t-1, 
                % representing the swap value at time t+1. The second term is the expected
                % present value of the coupon fraction cash-flow at time t+1 
                SwapTree(j-1)={(Probs*SwapTree{j}'+(FloatLeg-FixedLeg))'.*Bond};
                
                % adjust the backward induction procedure to determine the tree for a swaption portfolio
                swaption_op=1-(repmat((FwdDates==T(j-1)),1,Status(j-1)) & SwapTree{j-1}<=0);
                SwapTree(j-1)={SwapTree{j-1}.*swaption_op};

            end
            % cut the dead branches
            SwapTree(end)=[];

        premium = SwapTree{1}*1000;
        fun = sqrt(sum((bmk_swp_premium-premium).^2))/100;

    end

end