%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solutions to the Introduction to Matlab             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% To start with we clear the memory, so the Matlab workspace is empty
clear all
%
% Give MATLAB the path, where you want to work. 
% Otherwise, MATLAB doesn't know where to look for, and the only thing it knows are the 'built-in' functions ...
%
path(path,'your path goes here')
%
%% 2 Basic Commands

%a)

% Either you can use A = [1 3 2;5 0 2]; 
% or you specify each entry: A(1,1)= 1; ...
% or what I use below:

A= [1 3 2
5 0 2];

%b) 

% Matrix of zeros
V=zeros(4,2);

% Matrix of ones
W=ones(3,1);

% Identity matrix 
I= eye(3,3);

% Create a vector Z=[1,2,...100]'
Z=[1:1:100]';
Y=[1:4:100]';

% Saving the output of a MATLAB run:
%
% We want to save the output of our exercise in the directory P:\Lectures assisted\HS2013\Doing13\Session8, 
% in a fle called Session8.m
%
% We start by creating the following object:
DFILE=['your path goes here']; % ==> Output file
%and we save it to the P:\Lectures assisted\HS2013\Doing13\Session8 directory:

save(DFILE);
clear;

% You retrieve the results via the command load:
load 'your path goes here';

%c)
% An easy way to do this is:
% Column vector with 10 normally distributed random variables
% RAND([M,N]) returns an M-by-N matrix.

R0 = randn(10,1);

% Generate values from the uniform distribution on the interval [a, b].
%          r = a + (b-a).*rand(100,1);

R10 = 10+randn(1,10);

% Uniform distribution

U = rand(10,1);

% A nicer way is to use: 
% needs statistics toolbox!
% R_0 = random('Normal',0,1,10,1); 
% R_10= 10+random('Normal',0,1,1,10);
% U_1 = random('Uniform',0,1,10,1);

% You can find the probability distributions on
% http://www.mathworks.ch/ch/help/stats/random.html

% d)
A(1,2) = 8;
% e) 
B = A(2,:);
% f)
C = [4 1 2; 1 3 3];
D = C(C<=3);

%% 3 Basic Manipulations
clear all %clears the workspace
clc %clears the command window

%a)
E=[4 1;1 3]
F=[1 1;2 3]
[n,m]=size(E);

E+F  
E-F

%b)
% Matrix multiplication (not element by element):
G=E*F 
% Multiplication element by element:
H=E.*F
% Division element by element:
X=E./F

%c)
K=inv(E)*F
K=E^-1*F    
K=E\F    
% Determinant of a matrix:
det(K)

%% 4 Basic plotting and statistics
clear all %clears the workspace
clc %clears the command window

%a)
%x = random('Uniform',0,100,1,100);
x = 100*rand(100,1);

%b)
plot(x);

%c)
x1 = x(1:50,1);
x2 = x(51:100,1);

figure(1);
subplot(211);
plot(x1);
subplot(212);
plot(x2);

%other possibility
figure(2);
subplot(121);
plot(x1);
subplot(122);
plot(x2);

%editing
figure(3);
plot(x,'r:');
title('100 RANDOM DRAWS');
xlabel('number of the draw');
ylabel('drawn number');
legend('random draws');
axis([0 100 0 100]);

%d)
mean_x =mean(x);
std_x  =std(x);
mean_x1=mean(x1);
std_x1 =std(x1);
mean_x2=mean(x2);
std_x2 =std(x2);

%repeating b)
mean_x_vector=zeros(100,1);       
mean_x_vector(1:100,1)=mean_x;
figure(4);
hold on   ;                        
plot(x,'b');
plot(mean_x_vector,'r--');
title('100 RANDOM DRAWS');
xlabel('number of the draw');
ylabel('drawn number');
legend('random draws','Location','SouthOutside')   ;  
axis([0 100 0 100]);
hold off;

% You can find everything for editing plots on
% http://www.mathworks.ch/ch/help/matlab/creating_plots/using-high-level-plotting-functions.html

%% 5 Some easy loops
%a)
g=[0:0.1:10]'
size(g) %if you are unsure about the size of g
for j=1:101;          
    g(j)=g(j)+0.05;
end
g

%b)
A=2
while A<10;
    A=A+1
end

%c) IF-loop
r=rand(1,1);
if r<0.9;
    disp('This number is smaller than 0.9');
else 
    disp('This number is greater than 0.9');
end

%% 6 Geometric series
% a)
s=zeros(10,1);
q=0.5;
for n=1:10
    s(n,1)=q^(n-1);
end

%b)
SUM=sum(s);
disp('Sum for n=10, q=0.5')
disp(SUM)

%c)
s=zeros(10,3);
q1=0.2;
q2=0.5;
q3=0.8;
for n=1:10
    s(n,1)=q1^(n-1);
    s(n,2)=q2^(n-1);
    s(n,3)=q3^(n-1);
end
SUM=zeros(1,3);
for i=1:3
    SUM(1,i)=sum(s(:,i));
end
disp('Sums for q1=0.1, q2=0.5, q3=0.8')
disp(SUM)

%% 7 Difference equations
clear;

x=20;
T=(1:x)';
f=0.7;
y=zeros(x,1);
w=zeros(x,1);
y(1,1)=1.5;

% a)
% y_t = f*y_(t-1) + w_t
for t=1:x-1
    y(t+1,1)=f*y(t,1)+w(t+1,1);
end
figure;
plot(T,y,'b-');
% set f=1:
f=1;
y1=zeros(x,1);
y1(1,1)=1.5;
for t=1:x-1
    y1(t+1,1)=f*y1(t,1)+w(t+1,1);
end
figure;
plot(T,y1,'r-');

%b) Immediate and future effects: IRF
f=0.7;

% First approach: Analytical
IRF_analytical=zeros(x,1);
for t=1:x-1
    IRF_analytical(t+1,1)=f^(t-1);
end
figure
plot(T,IRF_analytical,'g:');

%Second Approach: Difference between shocked and unshocked path of y
w(2,1)=1;
y_unshocked=y;
y_shocked=zeros(x,1);
y_shocked(1,1)=1.5;
for t=1:x-1
    y_shocked(t+1,1)=f*y_shocked(t,1)+w(t+1,1);
end
IRF_difference=y_shocked-y_unshocked;
figure
plot(T,IRF_difference,'m:');

%c) 
% First approach: Analytical
IRF1_analytical=zeros(x,1);
for t=1:x-1
    IRF1_analytical(t+1,1)=sum(IRF_analytical(1:t+1,1));
end
figure
plot(T,IRF1_analytical,'g--');

%Second Approach: Difference between shocked and unshocked path of y
w(2:x,1)=1;
y_shocked1=zeros(x,1);
y_shocked1(1,1)=1.5;
for t=1:x-1
    y_shocked1(t+1,1)=f*y_shocked1(t,1)+w(t+1,1);
end
IRF_difference1=y_shocked1-y_unshocked;
figure
plot(T,IRF_difference1,'m--');

%d) Optional
w(2:x,1)=randn(x-1,1)*0.1;
ys=y;
for t=1:x-1
    ys(t+1,1)=f*ys(t,1)+w(t+1,1);
end
figure;
plot(T,ys,'b-');

%% 8 An almost infinite geometric series
clear all;
clc

sold=1;
q=0.5;
diff=1;
for n=1:1000000000; 
    snew=sold+q^n;
    diff=snew-sold; 
    if diff<=0.000000000001; 
        break;
    else sold=snew;
    end; 
end; 
n

%% 9 Long term effect
clear all;
clc

f=0.7;
yold=0;
for t=1:1000000000
    ynew=f^(t-1)+yold;
    diff=ynew-yold;
    if diff<=0.0000001;
        break;
    else yold=ynew;
    end;
end
t
disp('Long-Term Effect (after t periods)')
disp(ynew)
% NOT NECESSARY ACCORDING TO THE QUESTION
disp('Theoretical Value')
disp(1/(1-f))

%% 10 Function file

% Hypotenuse: Pythagoras sqrt(a^2+b^2)=c
a=input('Adjacent leg: ');                % German: Ankathete
b=input('Opposite leg: ');                % German: Gegenkathete
c=hypotenuse(a,b);                        % Call function-file "hypotenuse.m"
disp('Hypotenuse = Square Root of the Sum of the Adjacent- and the Opposite Leg')
disp(c)

%b)
disp('Press any key to draw a random number N(0,1)')
pause
z=randn(1,1);
judgement=draw(z);                      % Call function file "draw.m"
disp(judgement)
disp(z)