classdef Matrices < handle
    %Main matrices associated with the IOtables object model 
    %   
    %
    
    properties
                
    end
    
    properties ( SetAccess = private , Hidden = false )
    
       Parent
       Make
       Use
       FD
       Ab
       Ab_Dom
       Af
       Af_Dom
       Bb_Dom
       Bb
       Bf_Dom
       Bf
       FDShares
       Initialb
       Initialb_Dom
       Initialf
       Initialf_Dom
       InitialFD
       TotVA
       TotM
       TotFD
       TotIntermediates
          
    end
    
    properties (Dependent = true, SetAccess = private)
       
        
    end
    
    methods
        
        function obj = Matrices( MyParent )
            % Constructor that tracks what type of parent object has
            obj.Parent = MyParent;            
                        
            % Loading the matrices
            Year_ = obj.Parent.Year;
            %
            % - Make matrix                        
            Fname = strcat(obj.Parent.Settings.IOTablesPath,'\IONom\',obj.Parent.Settings.MakeName);
            Fname = strcat(Fname,'\',obj.Parent.Settings.MakeName,'_',int2str(Year_),'.dat');
            fid = fopen(Fname);       
            Make_ = fscanf( fid , '%f' , [obj.Parent.Settings.NoIndustries,obj.Parent.Settings.NoIndustries] );
            Make_ = Make_';
            obj.Make = Make_;
            fid = fclose(fid);
            %
            % - Use matrix            
            Fname = strcat(obj.Parent.Settings.IOTablesPath,'\IONom\',obj.Parent.Settings.UseName);
            Fname = strcat(Fname,'\',obj.Parent.Settings.UseName,'_',int2str(Year_),'.dat');
            fid = fopen(Fname);       
            Use_ = fscanf( fid , '%f' , [obj.Parent.Settings.NoIndustries+1,obj.Parent.Settings.NoIndustries+1] );
            Use_ = Use_';
            obj.Use = Use_;
            fid = fclose(fid);
            %
            % - FD matrix            
            Fname = strcat(obj.Parent.Settings.IOTablesPath,'\IONom\',obj.Parent.Settings.FDName);
            Fname = strcat(Fname,'\',obj.Parent.Settings.FDName,'_',int2str(Year_),'.dat');
            fid = fopen(Fname);       
            FD_ = fscanf( fid , '%f' , [obj.Parent.Settings.NoFD,obj.Parent.Settings.NoIndustries] );
            FD_ = FD_';
            obj.FD = FD_;
            fid = fclose(fid);                        
            
            % Getting rid of inventories, scrap, and second-hands
            FD_(:,134:141) = zeros(obj.Parent.Settings.NoIndustries,8);
            FD_(192:193,:) = zeros(2,size(FD_,2));
            Make_(192:193,:) = zeros(2,size(Make_,2));
            Make_(:,192:193) = zeros(size(Make_,1),2);
            Use_(192:193,:) = zeros(2,size(Use_,2));
            Use_(:,192:193) = zeros(size(Use_,1),2);
            
            %
            % - Distribute the markups to calculate the FDShare matrix
            FDShare =  FD_ ;
            % Markups on consumption expenditures
            FDShare(:,1:75) = DistMarkups( FDShare(:,1:75) , FDShare(:,76:82) );
            FDShare(:,76:82) = zeros(size(FDShare,1),7);
            % Markups on investment expenditures on equipment and software
            FDShare(:,84:110) = DistMarkups( FDShare(:,84:110) , FDShare(:,111:117) );
            FDShare(:,111:117) = zeros(size(FDShare,1),7);
            % Markups on investment in nonresidential structures
            FDShare(:,118) = DistMarkups( FDShare(:,118) , FDShare(:,119:125) );
            FDShare(:,119:125) = zeros(size(FDShare,1),7);
            % Markups on investment in residential structures
            FDShare(:,126) = DistMarkups( FDShare(:,126) , FDShare(:,127:133) );
            FDShare(:,127:133) = zeros(size(FDShare,1),7);
            % Markups on exports of goods
            FDShare(:,142) = DistMarkups( FDShare(:,142) , FDShare(:,143:149) );
            FDShare(:,143:149) = zeros(size(FDShare,1),7); 
            % Markups on exports of services
            FDShare(:,150) = DistMarkups( FDShare(:,150) , FDShare(:,151:157) );
            FDShare(:,151:157) = zeros(size(FDShare,1),7); 
            % Markups on imports of goods
            FDShare(:,158) = DistMarkups( FDShare(:,158) , FDShare(:,159:162) );
            FDShare(:,159:162) = zeros(size(FDShare,1),4);
            % Markups on Federal Government defense spending
            FDShare(:,164:165) = DistMarkups( FDShare(:,164:165) , FDShare(:,166:172) );
            FDShare(:,166:172) = zeros(size(FDShare,1),7);
            % Markups on Federal Government nondefense spending
            FDShare(:,173:174) = DistMarkups( FDShare(:,173:174) , FDShare(:,175:181) );
            FDShare(:,175:181) = zeros(size(FDShare,1),7);
            % Markups on State and Local government spending
            FDShare(:,182:183) = DistMarkups( FDShare(:,182:183) , FDShare(:,184:190) );
            FDShare(:,184:190) = zeros(size(FDShare,1),7);

            % Calculate total imports, final demand, intermediate demand,
            % and value added
            TotM = abs( FDShare(:,158) ) + abs(FDShare(:,163) );
            TotFD = sum( FDShare , 2 ) + TotM ;
            TotIntermediates = sum( Use_(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries) , 2 );
            TotVA = max( Use_(obj.Parent.Settings.NoIndustries+1,1:obj.Parent.Settings.NoIndustries) , 0 );
            TotVA = TotVA';
            
            obj.TotVA = TotVA;
            obj.TotFD = TotFD;
            obj.TotIntermediates = TotIntermediates;
            obj.TotM = TotM;
            
            % In each of these rows we find the fraction of output of an
            % industry that flows to a commodity (in column)
            TIndToComm = Make_;
            TIndToComm = kron(ones(1,size(TIndToComm,2)),(sum(TIndToComm,2)~=0))...
                         .*TIndToComm./kron(ones(1,size(TIndToComm,2)),(sum(TIndToComm,2)+(sum(TIndToComm,2)==0)));                     
            % In each of these rows we consider what fraction of commodity
            % output is used as intermediate input by industries
            TCommToInd = Use_(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries);
            TCommToInd = kron(ones(1,size(TCommToInd,2)),(sum([TotFD TotIntermediates],2)~=0))...
                         .*TCommToInd./kron(ones(1,size(TCommToInd,2)),(sum([TotFD TotIntermediates],2)+(sum([TotFD TotIntermediates],2)==0)));
            % In this matrix, we consider what fraction of commodity output
            % is used as final demand
            TCommToFD  = ( sum([TotFD TotIntermediates],2) ~= 0 ).* ...
                            TotFD./(sum([TotFD TotIntermediates],2)+(sum([TotFD TotIntermediates],2)==0));
            TCommToFD  = eye( obj.Parent.Settings.NoIndustries ).* ...
                    kron(ones(1,size(TCommToFD,1)),TCommToFD);
            % In this matrix, we consider what fraction of industry output
            % is used as intermediate inputs by another industry
            TIndToInd = TIndToComm*TCommToInd;
            % In this matrix, we consider what fraction of industry output
            % flows to final demand
            TIndToFD = TIndToComm*TCommToFD;
            % In this matrix, we consider what fraction of commodity
            % imports are used as intermediate inputs
            TMToInd = TCommToInd;
            % In this matrix, we consider what fraction of commodity
            % imports flow directly to final demand
            TMToFD  = TCommToFD;
        
            % -----------------------------------------------------------
            % This completes the matrices we need to construct the forward
            % Markov transition matrices and initial vectors. We construct
            % these here
            % Af moves [Y,M,FD]  to FD
            % Af_Dom moves [Y,FD] to FD
            % ------------------------------------------------------------
            obj.Af = zeros(3.*obj.Parent.Settings.NoIndustries,3.*obj.Parent.Settings.NoIndustries);
            obj.Af(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries) = TIndToInd';
            obj.Af(1:obj.Parent.Settings.NoIndustries,(1+obj.Parent.Settings.NoIndustries):(2.*obj.Parent.Settings.NoIndustries)) = TMToInd';
            obj.Af((1+2.*obj.Parent.Settings.NoIndustries):(3.*obj.Parent.Settings.NoIndustries),1:obj.Parent.Settings.NoIndustries) = TIndToFD';
            obj.Af((1+2.*obj.Parent.Settings.NoIndustries):(3.*obj.Parent.Settings.NoIndustries),(1+obj.Parent.Settings.NoIndustries):(2.*obj.Parent.Settings.NoIndustries)) = TMToFD'; 
            obj.Af_Dom = zeros(2.*obj.Parent.Settings.NoIndustries,2.*obj.Parent.Settings.NoIndustries);
            obj.Af_Dom(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries) = TIndToInd';            
            obj.Af_Dom((1+obj.Parent.Settings.NoIndustries):(2.*obj.Parent.Settings.NoIndustries),1:obj.Parent.Settings.NoIndustries) = TIndToFD';            
            
            % This is the matrix with the share in output of commodity i
            % that industry j makes
            TCommFromInd = Make_';
            TCommFromInd = kron(ones(1,size(TCommFromInd,2)),(sum(TCommFromInd,2)~=0))...
                         .*TCommFromInd./kron(ones(1,size(TCommFromInd,2)),(sum(TCommFromInd,2)+(sum(TCommFromInd,2)==0)));                     
            Dum = ( ones(size(TotM,1),1) - (TotM./(TotFD+TotIntermediates)) );
            Dum( isnan(Dum) ) = 0 ;
            TCommFromInd = TCommFromInd.*kron(ones(1,size(TCommFromInd,2)),Dum);
           
            TCommFromInd_Dom = Make_';
            TCommFromInd_Dom = kron(ones(1,size(TCommFromInd_Dom,2)),(sum(TCommFromInd_Dom,2)~=0))...
                         .*TCommFromInd_Dom./kron(ones(1,size(TCommFromInd_Dom,2)),(sum(TCommFromInd_Dom,2)+(sum(TCommFromInd_Dom,2)==0)));
            
            % This is row-vector with in each element the value added share
            % of each industry in gross output
            VAShare = TotVA'./(sum(Use_(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries),1)+TotVA');
            VAShare(isnan(VAShare)) = 0;
            
            % This is a matrix with (i,j)-th element being the intermediate
            % share in commodity i of industry j
            IntShare = Use_(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries)./ ...
                       kron(ones(obj.Parent.Settings.NoIndustries,1),(sum(Use_(1:obj.Parent.Settings.NoIndustries,1:obj.Parent.Settings.NoIndustries),1)+TotVA'));
            IntShare(isnan(IntShare)) = 0;
            
            % In this matrix we track what fraction of final demand comes
            % from value added of different industries
            TFDFromVA = TCommFromInd.*kron(ones(size(TCommFromInd,1),1),VAShare);
            TFDFromVA_Dom = TCommFromInd_Dom.*kron(ones(size(TCommFromInd_Dom,1),1),VAShare);
            
            % In this matrix we track what fraction of final demand comes
            % from intermediate inputs
            TFDFromInt = TCommFromInd*IntShare';
            TFDFromInt_Dom = TCommFromInd_Dom*IntShare';
            
            % In this matrix, we track what fraction of final demand comes
            % from imports.
            TFDFromM = eye(size(TotM,1)).* ...
                    kron(ones(1,size(TotM,1)),(TotM./(TotFD+TotIntermediates)));
            TFDFromM(isnan(TFDFromM)) = 0;
            
            % In this matrix, we track what fraction of intermediate inputs
            % (in commodities) comes from VA
            TIntFromVA = TFDFromVA ;
            TIntFromVA_Dom = TFDFromVA_Dom ;
            
            % In this matrix, we track what fraction of intermediate inputs
            % (in commodities) comes from intermediate inputs
            TIntFromInt = TFDFromInt;
            TIntFromInt_Dom = TFDFromInt_Dom;
            
            % In this matrix, we trach what fraction of intermediate inputs
            % (in commodities) comes from imports 
            TIntFromM = TFDFromM;
            
            NoInd = obj.Parent.Settings.NoIndustries;
                                                
            % -----------------------------------------------------------
            % This completes the matrices we need to construct the backward
            % Markov transition matrices and initial vectors. We construct
            % these here
            % Ab moves [FD,Y,VA,M] to M and VA
            % Ab_Dom moves [ FD , Y , VA ] to VA
            % -----------------------------------------------------------
            obj.Ab = zeros(4.*NoInd,4.*NoInd);
            obj.Ab((NoInd+1):(2.*NoInd),1:NoInd) = TFDFromInt';
            obj.Ab((2.*NoInd+1):(3.*NoInd),1:NoInd) = TFDFromVA';
            obj.Ab((3.*NoInd+1):(4.*NoInd),1:NoInd) = TFDFromM';
            obj.Ab((NoInd+1):(2.*NoInd),(NoInd+1):(2.*NoInd)) = TIntFromInt';
            obj.Ab((2.*NoInd+1):(3.*NoInd),(NoInd+1):(2.*NoInd)) = TIntFromVA';
            obj.Ab((3.*NoInd+1):(4.*NoInd),(NoInd+1):(2.*NoInd)) = TIntFromM';
            
            obj.Ab_Dom = zeros(3.*NoInd,3.*NoInd);
            obj.Ab_Dom((NoInd+1):(2.*NoInd),1:NoInd) = TFDFromInt_Dom';
            obj.Ab_Dom((2.*NoInd+1):(3.*NoInd),1:NoInd) = TFDFromVA_Dom';
            obj.Ab_Dom((NoInd+1):(2.*NoInd),(NoInd+1):(2.*NoInd)) = TIntFromInt_Dom';
            obj.Ab_Dom((2.*NoInd+1):(3.*NoInd),(NoInd+1):(2.*NoInd)) = TIntFromVA_Dom';
              
            B = FDShare;
            B(: ,[ 158 163 ]) = zeros( size(B,1) , 2 );
            
            obj.Bb = B./kron(sum(B,1),ones(size(B,1),1));
            obj.Bb(isnan(obj.Bb))=0;
            obj.Bb_Dom = [ obj.Bb ; zeros(2.*NoInd,size(obj.Bb,2)) ];
            obj.Bb = [ obj.Bb ; zeros(3.*NoInd,size(obj.Bb,2)) ];
            
            obj.Bf = B./kron(sum(B,2),ones(1,size(B,2)));
            obj.Bf(isnan(obj.Bf))=0;
            obj.InitialFD = obj.Bf'*TotFD;
            obj.Bf_Dom = [ obj.Bf' zeros(size(obj.Bf,2),NoInd) ];
            obj.Bf = [ obj.Bf' zeros(size(obj.Bf,2),2.*NoInd) ];
                        
            obj.Initialb = [ TotFD ; zeros(3.*NoInd,1) ];
            obj.Initialb_Dom = [ TotFD ; zeros( 3.*NoInd,1) ];
            obj.Initialf = [ TotVA ; TotM ; zeros(NoInd,1) ];
            obj.Initialf_Dom = [ (TotVA+TotM) ; zeros(NoInd,1) ];
            
        end
        
        % Methods used to change the parameter values
        
    end
    
end

