

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%  Auxiliary Matlab codes %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('final_raw.mat')
% Apply the one-sided HP filter
lambda = 1600; % smoothing parameter, may need to tune it up or down 
% but we use 1600 for quarterly series
trend = hpfilter(gdp, lambda);

% Plot the original data and the trend component
figure;
plot(year, gdp, 'b', 'LineWidth', 2);
hold on;
plot(year, trend, 'r', 'LineWidth', 2);
title('Original Data and One-Sided HP Filtered Trend');
xlabel('Year');
ylabel('GDP(Billions of NOK)');
legend('Original Data', 'One-Sided HP Filtered Trend');
hold off;

%%%%%%%%%%%%%%%%
One-sided hpfilter
%%%%%%%%%%%%%%%%%

%% GDP

lambda=1600; %% Being quarterly series

% One-sided HP filter
[TTbl1S_gdp, CTbl1S_gdp, h1S] = hpfilter(gdp_olvar, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(year, gdp_olvar, 'b', 'LineWidth', 2);
hold on;
plot(year, TTbl1S_gdp, 'r', 'LineWidth', 2);
hold on;
plot(year, CTbl1S_gdp, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('GDP (Millions of NOK at 2005 Prices)');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;

%% Emission

lambda=1600 %% Being quarterly series

% One-sided HP filter
[TTbl1S_emiss, CTbl1S_emiss, h1S] = hpfilter(emiss, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(year, emiss, 'b', 'LineWidth', 2);
hold on;
plot(year, TTbl1S_emiss, 'r', 'LineWidth', 2);
hold on;
plot(year, CTbl1S_emiss, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('Emission (1000 Tonnes CO2 Eqv)');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;

%% Rate of Returns on Long-Term Govt Bonds
lambda=1600 %% Being quarterly series

% One-sided HP filter
[TTbl1S_bond, CTbl1S_bond, h1S] = hpfilter(bond, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(year, bond, 'b', 'LineWidth', 2);
hold on;
plot(year, TTbl1S_bond, 'r', 'LineWidth', 2);
hold on;
plot(year, CTbl1S_bond, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('(%) Rate of Returns on Govt. Bonds (10 Yrs)');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;

%% CPI
lambda=1600 %% Being quarterly series

% One-sided HP filter
[TTbl1S_cpi, CTbl1S_cpi, h1S] = hpfilter(cpi, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(year, cpi, 'b', 'LineWidth', 2);
hold on;
plot(year, TTbl1S_cpi, 'r', 'LineWidth', 2);
hold on;
plot(year, CTbl1S_cpi, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('CPI (Base 2015 = 100');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;


%%%%%%%%%%%%%%%%
kalman_filter
%%%%%%%%%%%%%%%%%
load('final_raw.mat')
% Define the state space model
A = 1;  % State transition matrix (assumed identity matrix for simplicity)
H = 1;  % Measurement matrix
Q = 0.1;  % Process noise covariance
R = 1;  % Measurement noise covariance

% Initial state estimate and covariance
x_hat = 0;  % Initial state estimate (you can set this based on your prior knowledge)
P = 1;  % Initial state covariance

% Generate synthetic data
true_gdp = gdp;  % True GDP (you can replace this with your actual data)
observed_gdp = true_gdp + 0.5 * randn(1, 33);  % Observed GDP with noise

% Kalman filter initialization
num_samples = length(observed_gdp);
filtered_gdp = zeros(1, num_samples);
for k = 1:num_samples
    % Prediction step
    x_hat_minus = A * x_hat;
    P_minus = A * P * A' + Q;

    % Update step
    K = P_minus * H' / (H * P_minus * H' + R);
    x_hat = x_hat_minus + K * (observed_gdp(k) - H * x_hat_minus);
    P = (1 - K * H) * P_minus;

    % Save the filtered estimate
    filtered_gdp(k) = x_hat;
end

% Plot the results
figure;
plot(1:num_samples, true_gdp, 'b-', 'LineWidth', 2, 'DisplayName', 'True GDP');
hold on;
plot(1:num_samples, observed_gdp, 'r.', 'DisplayName', 'Observed GDP');
plot(1:num_samples, filtered_gdp, 'g-', 'LineWidth', 2, 'DisplayName', 'Filtered GDP');
legend('Location', 'Best');
title('Kalman Filter for GDP Estimation');
xlabel('Time');
ylabel('GDP');
grid on;
hold off;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Finding demeaned series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example yearly data
yearlyData = rand(5, 1); % replace this with actual file or variable
% Take the logarithm of the yearly data
logYearlyData = log(yearlyData);

% Demean the log-transformed data
demeanedLogYearlyData = logYearlyData - mean(logYearlyData);

% Display original, log-transformed, and demeaned data
disp('Yearly Data:')
disp(yearlyData)

disp('Log-Transformed Yearly Data:')
disp(logYearlyData)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Breakpoint (BEAST Analysis)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('finalDataRaw.mat')
load('year.mat') %% works only with 'datetime' not with'timetable'

%%gdp
out = beast(GDP, 'start', 1990, 'deltat', 1/4)  % specifying start time and indicating that it is quarterly 
 printbeast(out)
 plotbeast(out)

 %%cpi
out = beast(CPI, 'start', 1990, 'deltat', 1/4)  % specifying start time and indicating that it is quarterly 
 printbeast(out)
 plotbeast(out)
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Filterd Variables vis a vis Historical series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% must load 'variables' plus 'year' before execution after estimation
%% can't close estimation results if not saved before running the following codes
figure;
plot(year, Y, 'b', 'LineWidth', 2); %% can replace E with TTbl1S_emiss from hp_filtered trend
hold on;
plot(year, oo_.FilteredVariables.E, 'r', 'LineWidth', 2);
title('Original series and Smoothed series');
xlabel('Year');
ylabel('GDP');
legend('Original Series', 'Filtered Series');
hold off;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Smoothed Variables vis a vis Historical series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
plot(year, E, 'b', 'LineWidth', 2); %% E=emissions
hold on;
plot(year, oo_.SmoothedVariables.E, 'r', 'LineWidth', 2);
title('Original series and Smoothed series');
xlabel('Year');
ylabel('Emissions');
legend('Original Series', 'Filtered Series');
hold off;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Descriptive statistics of Empirical time series
%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Compute descriptive statistics
statistics = [mean(gdp), std(gdp), skewness(gdp), kurtosis(gdp), min(gdp), max(gdp);
              mean(emiss), std(emiss), skewness(emiss), kurtosis(emiss), min(emiss), max(emiss);
              mean(cpi), std(cpi), skewness(cpi), kurtosis(cpi), min(cpi), max(cpi);
              mean(bond), std(bond), skewness(bond), kurtosis(bond), min(bond), max(bond)];

% Create a table to display the results
variable_names = {'Real GDP', 'Emissions', 'CPI', 'Real Interest Rate'};
stat_names = {'Mean', 'Std', 'Skewness', 'Kurtosis', 'Min', 'Max'};
result_table = array2table(statistics, 'VariableNames', stat_names, 'RowNames', variable_names);

% Display the table
disp(result_table);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Bad plots the following gave, but no errors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load the .mat file
load('final_raw.mat');

% Assuming 'year' is a column vector and 'gdp' is a column vector

% Generate quarterly 'year' values
quarters = linspace(year(1), year(end), numel(year)*4).';

% Interpolate quarterly GDP values
gdp_quarterly = interp1(year, gdp, quarters);

% Save the new variables to a new .mat file
save('your_quarterly_data_file.mat', 'quarters', 'gdp_quarterly');
plot(quarters,gdp_quarterly)





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
import from excel From practice folder, under Data Analysis, in Thesis
%%%%%%%%%%%%%%%%%%%%%%%%%%%

tQ=readtable('09190_20231209-112708_gdp.xlsx');
tQ.year=datetime(tQ.year,'InputFormat','uuuuQQQ','Format','uuuuQQQ');
head(tQ)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%importing from csv => same as in excel
%%%%%%%%%%%%%%%%%%%%%%%

date = tQ.DATE; % replace this with your datetime data
values = tQ.CPI; % replace this with your double data

% Plotting
figure;
plot(date, values, 'o-');
title('CPI in Norway');
xlabel('Year');
ylabel('CPI');
datetick('x', 'yyyyQQ', 'keepticks'); % Format the x-axis as per your preference
grid on;


%%%%%%%%%%%%%%%%%%%
%%transform double to datetime
%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Emissions

tr=readtable('13932_20231209-174139.csv');
tr.Properties.VariableNames(4) = "Emission";
tr = removevars(tr, ["industry","pollutant"]);

% Assuming your data is stored in a matrix (as table) named 'tR' and year column is
% present
tr.year = datetime(tr.year, 1, 1);


%%% after that,

% Assuming 'year' is the datetime column and 'Emission' is the double column
% Convert the data to timetable
tt = table2timetable(tr, 'RowTimes', 'year');

% Create a quarterly time vector from the start of 1990 to the end of 2022
quarterlyTime = datetime(1990, 1, 1):calquarters(1):datetime(2022, 12, 31);

% Expand the Emission values to match the quarterly time vector
expandedEmission = repelem(tt.Emission, 4);

% Create a new timetable with quarterly time vector and expanded Emission values
ttQuarterly = timetable(quarterlyTime', expandedEmission, 'VariableNames', {'Emission'});

% Display the resulting timetable
disp(ttQuarterly);

%%% and finally
ttQuarterly.Time=datetime(ttQuarterly.Time,'InputFormat','uuuuQQQ','Format','uuuuQQQ');

%% to change timetable to table
T = timetable2table(ttQuarterly)
%%or


%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Final
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('fredgraph.mat') %% I collected all the 4 series into one .mat file
fredgraph.Properties.VariableNames(1) = "year"; %% changed period to year
fredgraph.year=datetime(fredgraph.year,'InputFormat','uuuuQQQ','Format','uuuuQQQ'); %% changed to quarters




%%%%%%%%%%%%%%%%
% One-sided hpfilter
%%%%%%%%%%%%%%%%%

%% GDP

lambda=1600; %% Being quarterly series

% One-sided HP filter
[TTbl1S_gdp, CTbl1S_gdp, h1S] = hpfilter(GDP, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(fredgraph.year, GDP, 'b', 'LineWidth', 2);
hold on;
plot(fredgraph.year, TTbl1S_gdp, 'r', 'LineWidth', 2);
hold on;
plot(fredgraph.year, CTbl1S_gdp, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('GDP (Millions of Chained 2010 Euros)');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;

%% Emission

lambda = 1600;

% One-sided HP filter

[TTbl1S_emiss, CTbl1S_emiss, h1S] = hpfilter(Emission, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(fredgraph.year, Emission, 'b', 'LineWidth', 2);
hold on;
plot(fredgraph.year, TTbl1S_emiss, 'r', 'LineWidth', 2);
hold on;
plot(fredgraph.year, CTbl1S_emiss, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('Emission (1000 Tonnes CO2 Eqv)');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;

%% Rate of Returns on Long-Term Govt Bonds

% One-sided HP filter
[TTbl1S_bond, CTbl1S_bond, h1S] = hpfilter(BondYield, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(fredgraph.year, BondYield, 'b', 'LineWidth', 2);
hold on;
plot(fredgraph.year, TTbl1S_bond, 'r', 'LineWidth', 2);
hold on;
plot(fredgraph.year, CTbl1S_bond, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('(%) Rate of Returns on Govt. Bonds (10 Yrs)');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;

%% CPI

% One-sided HP filter
[TTbl1S_cpi, CTbl1S_cpi, h1S] = hpfilter(CPI, 'FilterType', 'one-sided', 'Smoothing', lambda);

% Plot the original data and the trend component
figure;
plot(fredgraph.year, CPI, 'b', 'LineWidth', 2);
hold on;
plot(fredgraph.year, TTbl1S_cpi, 'r', 'LineWidth', 2);
hold on;
plot(fredgraph.year, CTbl1S_cpi, 'g', 'LineWidth', 2);

title('Original Series and One-Sided HP Filtered Series');
xlabel('Year');
ylabel('CPI (Base 2015 = 100');
legend('Original Data', 'One-Sided HP Filtered Trend', 'One-Sided HP Filtered Cycles');
hold off;
&&&&&&&&&

%%%%%%%%%%%%%%
%%for filtered emission trend
%%%%%%%%%%%%%%%%%%
figure;
plot(fredgraph.year, TTbl1S_emiss, 'b', 'LineWidth', 2);

%%%%%%%%%%%%%%
%%for filtered emission cycle
%%%%%%%%%%%%%%%%%%
figure;
plot(fredgraph.year, CTbl1S_emiss, 'b', 'LineWidth', 2);

%%%%%%%%%%%
%%Final data plot (deviation from steady states
%%%%%%%%%%%%%%%%%%%%%%%%%

%%GDP and Emission
figure;
plot(fredgraph.year, Y, 'b', 'LineWidth', 2);
hold on;
plot(fredgraph.year, E, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('GDP(Billions of NOK)');
legend('DGP', 'Emission');
hold off;

%% CPI
plot(fredgraph.year, PI, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('% Point Change in Inflation');
legend('Inflation');

%% RB

plot(fredgraph.year, RB, 'g', 'LineWidth', 2);
xlabel('Year');
ylabel('% Point Change in Interest Rate');
legend('Interest Rate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% for all figs. saving as .png


FolderName = ('G:\Image proccessing\test');   % using my directory
FigList = findobj(allchild(0), 'flat', 'Type', 'figure');
for iFig = 1:length(FigList)
  FigHandle = FigList(iFig);
  FigName   = num2str(get(FigHandle, 'Number'));
  set(0, 'CurrentFigure', FigHandle);
%   saveas(FigHandle, strcat(FigName, '.png'));
  saveas(FigHandle, fullfile(FolderName,strcat(FigName, '.png'))); % specify the full path
end


%another method%

% Specify the folder to save the figures
saveFolder = '/Users/binod/Desktop/thesis/graphs handdrawn/dataaa/practice/Final mode folder/figures';

% Create the folder if it doesn't exist
if ~exist(saveFolder, 'dir')
    mkdir(saveFolder);
end

% Save all open figures in the specified formats
figHandles = findobj('Type', 'figure');

for i = 1:numel(figHandles)
    % Save as JPG
    saveas(figHandles(i), fullfile(saveFolder, ['figure_', num2str(i), '.jpg']));
    
    % Save as FIG
    saveas(figHandles(i), fullfile(saveFolder, ['figure_', num2str(i), '.fig']));
    
    % Save as PNG
    saveas(figHandles(i), fullfile(saveFolder, ['figure_', num2str(i), '.png']));
end
%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%
%% To reproduce shock decomposition

% Access the shock decomposition results
shock_decomposition = oo_.shock_decomposition;

% Specify the endogenous variable index and time period
endogenous_variable_index = 1;  % Change this to the index of the desired endogenous variable
time_period = 1:120;            % Change this to the desired time period

% Extract the contributions of shocks and initial conditions for the specified variable and time period
contributions = squeeze(shock_decomposition(endogenous_variable_index, 1:M_.exo_nbr+1, time_period));

% Get the names of shocks and initial conditions
shock_names = M_.exo_names;
initial_condition_name = 'Initial Condition';

% Create a bar plot
figure;
bar(time_period, contributions, 'stacked');
legend([shock_names; initial_condition_name], 'Location', 'bestoutside');
title(['Shock Decomposition for ', M_.endo_names{endogenous_variable_index}]);
xlabel('Time Period');
ylabel('Contribution');

%Note: remember to change 'datetime' to 'numeric double' before inserting year
%along x-axis

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check for NaN values in Y
missingValues = isnan(Y);

Y_fixed = Y;
Y_fixed(missingValues) = 0;
mean(Y_fixed)

% Replace NaN values with a specific value (e.g., 0)
PI_fixed = PI;
PI_fixed(missingValues) = 0;

% Replace NaN values with a specific value (e.g., 0)
E_fixed = E;
E_fixed(missingValues) = 0;

% Replace NaN values with a specific value (e.g., 0)
RB_fixed = RB;
RB_fixed(missingValues) = 0;
mean(RB_fixed) % observe mean


%%%%%%%%%%%%%%%%%%%%%%%%%%%
Some additional codes
%%%%%%%%%%%%%%%%%%%%%%%%%%

                  
figure;
plot(year, oo_.UpdatedVariables.Mean.PI, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.PI, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Inflation');
legend('Updated Inflation', 'Forecasted Inflation');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.RB, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.RB, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Returns on Bonds');
legend('Updated Returns on Bonds', 'Forecasted Returns on Bonds');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.E, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.E, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Emisssion');
legend('Updated Emisssion', 'Forecasted Emisssion');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.IP, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.IP, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Private Investment');
legend('Updated Private Investment', 'Forecasted Private Investment');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.IG, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.IG, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Public Investment');
legend('Updated Public Investment', 'Forecasted Public Investment');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.W, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.W, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Wage Rate');
legend('Updated Wage Rate', 'Forecasted Wage Rate');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.L, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.L, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Labor');
legend('Updated Labor', 'Forecasted Labor');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.R, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.R, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Returns to Capital');
legend('Updated Returns to Capital', 'Forecasted Returns to Capital');
hold off;

figure;
plot(year, oo_.UpdatedVariables.Mean.MC, 'b', 'LineWidth', 2);
hold on;
plot(year, oo_.FilteredVariables.Mean.MC, 'r', 'LineWidth', 2);
xlabel('Year');
ylabel('Marginal Cost');
legend('Updated Marginal Cost', 'Forecasted Marginal Cost');
hold off;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   End of coding
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


