clear;
clc;

% The BEA IO data is more disaggregated than the BEA-BLS Production Account
% Use the aggregated productivity data for the more disaggregated industry

% All output files sort the industries in the following order:
% 1 .Farms
% 2. Forestry, fishing, and related activities
% 3. Oil and gas extraction
% 4. Mining, except oil and gas
% 5. Support activities for mining
% 6. Utilities
% 7. Construction
% 8. Wood products
% 9. Nonmetallic mineral products
% 10. Primary metals
% 11. Fabricated metal products
% 12. Machinery
% 13. Computer and electronic products
% 14. Electrical equipment, appliances, and components
% 15. Motor vehicles, bodies and trailers, and parts
% 16. Other transportation equipment
% 17. Furniture and related products
% 18. Miscellaneous manufacturing
% 19. Food and beverage and tobacco products
% 20. Textile mills and textile product mills
% 21. Apparel and leather and allied products
% 22. Paper products
% 23. Printing and related support activities
% 24. Petroleum and coal products
% 25. Chemical products
% 26. Plastics and rubber products
% 27. Wholesale trade
% 28. Motor vehicle and parts dealers
% 29. Food and beverage stores
% 30. General merchandise stores
% 31. Other retail
% 32. Air transportation
% 33. Rail transportation
% 34. Water transportation
% 35. Truck transportation
% 36. Transit and ground passenger transportation
% 37. Pipeline transportation
% 38. Other transportation and support activities
% 39. Warehousing and storage
% 40. Publishing industries, except internet (includes software)
% 41. Motion picture and sound recording industries
% 42. Broadcasting and telecommunications
% 43. Data processing, internet publishing, and other information services
% 44. Federal Reserve banks, credit intermediation, and related activities
% 45. Securities, commodity contracts, and investments
% 46. Insurance carriers and related activities
% 47. Funds, trusts, and other financial vehicles
% 48. Housing
% 49. Other real estate
% 50. Rental and leasing services and lessors of intangible assets
% 51. Legal services
% 52. Computer systems design and related services
% 53. Miscellaneous professional, scientific, and technical services
% 54. Management of companies and enterprises
% 55. Administrative and support services
% 56. Waste management and remediation services
% 57. Educational services
% 58. Ambulatory health care services
% 59. Hospitals
% 60. Nursing and residential care facilities
% 61. Social assistance
% 62. Performing arts, spectator sports, museums, and related activities
% 63. Amusements, gambling, and recreation industries
% 64. Accommodation
% 65. Food services and drinking places
% 66. Other services, except government
% ***************************************************************************/

% Define the raw data directory path
rawDataDir = fullfile(cd, '..', '..', 'data_raw', 'productivity'); % Path to raw data directory

% Define file paths
xlsFilePath = fullfile(rawDataDir, 'BEA-BLS-industry-level-production-account-1987-2019.xlsx'); % Load from data_raw
matFilePath = fullfile(cd, 'MFP.mat'); % Save in current directory (data_clean/productivity)

% Don't import Excel files if the .mat file already exists in the current directory
if exist(matFilePath, 'file') ~= 2  
    count = 1;

    % Import the frequency of price changes across industries
    MFP = xlsread(xlsFilePath, 'Integrated MFP Index', 'B3:AH63');
    log_MFP = log(MFP);

    % Disaggregate "Retail trade" into finer industry categories
    log_MFP = productivity_disaggregate(log_MFP, 28, 4);

    % Disaggregate "Real estate" into "Housing" and "Other real estate"
    log_MFP = productivity_disaggregate(log_MFP, 48, 2);

    % Disaggregate "Hospitals and nursing and residential care facilities"
    log_MFP = productivity_disaggregate(log_MFP, 59, 2);

    % Interpolate annual data to quarterly frequency
    n = size(log_MFP,1);
    t = 1:size(log_MFP,2);
    log_MFP_quarterly = zeros(n, 4*(size(log_MFP,2)-1) + 1);

    for i = 1:n
        log_MFP_quarterly(i,:) = interp1(t, log_MFP(i,:), 1:0.25:length(t));
    end

    % Detrending log productivity sequences
    log_MFP_linear_detrended = detrend(log_MFP_quarterly');

    % Compute variance-covariance matrix of shocks
    Sigmaz = cov(log_MFP_linear_detrended);

    % Clear temporary variables to free memory
    clear count MFP log_MFP n t log_MFP_quarterly i

    % Save the processed data in the current directory (data_clean/productivity)
    save(matFilePath, 'log_MFP_linear_detrended', 'Sigmaz');

else
    % Load preprocessed data if it already exists in the current directory
    load(matFilePath);
end