[theta]=unimcmctest

function [theta]=unimcmctest
%Fits the univariate Normal Inverse Gaussian distribution using MCMC Gibbs
%sampler
%Author: Panagiotis Braimakis

%Number of returns
n=length(y);

%MCMC sampling scheme according to the GIBBS sampler and by using the
%function randraw
p=input('How many iterations do you want? ');
clc;
disp('Calculating');

%Space allocation
a=zeros(p,1);
b=zeros(p,1);
m=zeros(p,1);
d=zeros(p,1);
gg=zeros(p,1);

%%%%%%%%%%%%%%%%%%%%%SEEDS%%%%%%%%%%%%%%%%%%%
a(1)=5;   %for 
b(1)=5;   %for 
m(1)=0;   %for 
d(1)=1;   %for 
gg(1)=sqrt(a(1)^2-b(1)^2);

theta2=zeros(2,1); %the prior mean vector
%theta2=[2;1];
C2=[10 0;0 10]; % the prior covariance matrix
%Likelihood parameters
%the design matrix
X=[ones(n,1)];
Y=y;

%hyperparameters a0,a1,a2,a3,a4
a0=0;
a3=0.001;
a4=0.001;

%START the simulations (GIBBS sampler)
for i=1:p

%1. zi~GIG(-1,(q(yi))^1/2,a) ; i=1,...,n
    q=1+((y-m(i))./d(i)).^2; %for q(yi)
        
    for j=1:length(q)
    gig1=(d(i)*sqrt(q(j)))^2;
    gig2= a(i)^2;
    z(i,j)=randraw('gig',[-1 ,gig1,gig2],1);
    end

%2. Simulation from the posterior multivariate normal with d and D
%prior parameters               
A1=[X,z(i,:)'];
%the known covariance matrix size nxn
C1=diag(z(i,:));

%the data return matrix
dd=A1' * inv(C1) * Y + inv(C2) * theta2;
D=inv(A1' * inv(C1)* A1 + inv(C2));
bd=mvnrnd(D*dd,D,1);
m(i+1)=bd(1);
b(i+1)=bd(2);

%3 Gamma-TN scheme
%^2=Gamma(a+n/2,+n*zbar/2+1/(2*theta)(w^2-(w+n*theta)^2/(1+n*zbar*theta))

%Parameters
ga(1)=(a0+n+1)/2 ;
ga(2)=a4+n*mean(1./z(i,:))/2-((a0 +n)^2)/(4*a3+2*n*mean(z(i,:))) ;  

%One random draw for d^2 
dsq=gamrnd(ga(1), 1/ga(2));
d(i+1)=sqrt(dsq);

%|~TN((w+n*theta)/(1+n*zbar*theta)*d,theta/(1+n*zbar*theta)) truncated
%from 0 to inf
%Parameters
nt(1)=(a0*d(i+1)+n*d(i+1))/(2*a3+n*mean(z(i,:))) ;
nt(2)= 1/(2*a3+n*mean(z(i,:)));

%One random draw for g
gg(i)=randraw('normaltrunc',[0,1e5,nt(1),(nt(2))^1/2],1);
a(i+1)=sqrt(gg(i)^2+b(i+1)^2);

%Store the parameter vectors
theta(i,:)=[a(i+1) b(i+1) m(i+1) d(i+1)];
i
[i a(i+1) b(i+1) m(i+1) d(i+1)]
if mod(i,100)==0

    plot(log(sum(z,2)))
    pause(0.005);
end

save pb2.mat theta z



end