function [rate_all, rate_ig, rate_eg, rate_race, price] = simulate_adoption_w_specific_net_ext_gss(tie, nrig, alpha, gamma, delta, k, init_price, max_time)

% function [rate_all, rate_group, rate_race, price] = simulate_adoption_w_specific_net_ext_gss(tie, nrig, alpha, gamma, delta, k, init_price, max_time)
% Simulates internet adoption rates for a given network and reservation
% price parameters; takes into account specific network externalities where 
% the reservation price is a function of the adoption rates in an individual's
% network; returns adoption rates and prices over time
% INPUT:    tie         - A tie matrix defining a network
%           nrig        - Node #, Race, Income, Group Id Matrix
%           alpha       - Exponential constant for the network effect
%           gamma       - Exponential constant for the income effect
%           delta       - Multiplicative constant for the network effect
%           k           - Multiplicative constant for the income effect
%           init_price  - Initial monthly price of the internet
%           max_time    - Number of time iterations in the simulation
% 
% Reservation price is computed as follows:
% Res Price = Income Effect  + Network Effect
% Res Price = k*income^gamma + delta* income^gamma * rate_net(1,time-1)^alpha
% 
% OUTPUT:   rate_all    - Matrix of adoption rates over time
%           rate_ig     - Matrix of adoption rates by income group over time
%           rate_eg     - Matrix of adoption rates by education group over time
%           rate_race   - Matrix of adoption rates by race over time
%           price       - Vector of internet prices over time
%
% Paul DiMaggio and Filiz Garip (c)
% Department of Sociology
% Princeton University, Harvard University
% Created: 28 March 2005
% Last Updated: 10 June 2008, 2:30pm
% File Name: simulate_adoption_w_specific_net_ext_gss.m

% Access Price Parameters 

p_eqm           = 28.74;                % Equilibrium Price Level - Estimated using off-peak prices in OECD data
a               = 3.34;                 % Rate of price decline 

max_obs         = size(tie,1);
adopt_ind       = zeros(max_obs,max_time);                  % Adoption indicator matrix - Rows: Nodes 1-N - Cols: Time 1-T 
adopt_in_net    = zeros(max_obs,max_time);                  % Number of adopters from the network of each node, at each time period
rate_in_net     = zeros(max_obs,max_time);                  % Adoption rate from the network of each node, at each time period
adopt_all       = zeros(1, max_time);                       % Overall number of adopters at each time period
rate_all        = zeros(1, max_time);                       % Overall adoption rate at each time period

net_effect      = zeros(max_obs,1);                         % Initialize network effect w/ zeros
price           = init_price * ones(1,max_time);            % Initialize prices
res_price       = zeros(max_obs,1);                         % Initialize reservation prices

income          = nrig(:,3);                                          
inc_effect      = k * income.^gamma;                        % Income effect on reservation price, constant over time 
net_size        = sum(tie,2);                               % Compute network size for each node


% Add a random perturbation to the reservation price of each individual
% that is distributed as normal with mean 0, and std dev $12.5 (half the 
% std deviation of the reservation prices without the perturbation term 
% under network externalites). 

pert = 12.5.*randn(size(nrig,1),1);

for t = 2:max_time
    
    net_effect              = income .^ gamma .* delta .* rate_in_net(:,t-1) .^ alpha;  % Network effect on reservation price
    res_price               = net_effect + inc_effect + pert;                           % Reservation price = Income effect + Network effect + Random Perturbation
    price(1,t)              = price(1,t-1)+ (a * rate_all(1,t-1) * (p_eqm - price(1,t-1)))/12;     % Update price    
    adopt_ind(:,t)          = (res_price >= price(1,t));                                % Update the adopters
    
    for i=1:max_obs                                                 % For each node    
        [I J]               = find(tie(:,i) == 1);                  % Identify the personal network
        adopt_in_net(i,t)   = sum(adopt_ind(I,t));                  % Compute the number of adopters within the network 
        if net_size(i,1) ~= 0
            rate_in_net(i,t)    = adopt_in_net(i,t) / net_size(i,1);    % Compute the adoption rate within the network
        end
    end
    
    adopt_all(1,t)          = sum(adopt_ind(:,t));                  % Count overall number of adopters
    rate_all(1,t)           = adopt_all(1,t)/max_obs;               % Compute overall rate of adopters
end

% Compute the adoption rates by each income-education group

nig             = size(unique(nrig(:,5)),1);            % Count number of income groups
neg             = size(unique(nrig(:,6)),1);            % Count number of education groups

adopt_ig        = zeros(nig, max_time);                 % Create a number of adopters by income group matrix- Rows: Groups - Columns: Time
rate_ig         = zeros(nig, max_time);                 % Create an adoption rate by income group matrix
igs             = zeros(nig,1);                         % Create a vector of income group sizes

adopt_eg        = zeros(neg, max_time);                 % Create a number of adopters by education group matrix- Rows: Groups - Columns: Time
rate_eg         = zeros(neg, max_time);                 % Create an adoption rate by education group matrix
egs             = zeros(neg,1);                         % Create a vector of education group sizes

% Compute adoption rates by income group

for i=1:nig                                             % For each income group
    [I J]               = find(nrig(:,5)==i);           % Find the group members
    igs(i)              = size(I,1);                    % Count the number of group members
    adopt_ig(i,:)       = sum(adopt_ind(I,:));          % Count the number of adopters over time
    rate_ig(i,:)        = adopt_ig(i,:)./igs(i);        % Compute the rate of adopters over time
end

% Compute adoption rates by education group

for i=1:neg                                             % For each education group
    [I J]               = find(nrig(:,6)==i-1);         % Find the group members
    egs(i)              = size(I,1);                    % Count the number of group members
    adopt_eg(i,:)       = sum(adopt_ind(I,:));          % Count the number of adopters over time
    rate_eg(i,:)        = adopt_eg(i,:)./egs(i);        % Compute the rate of adopters over time
end

% Compute the adoption rates by race

adopt_race   = zeros(2, max_time);                      % Create a number of adopters by race
rate_race    = zeros(2, max_time);                      % Create an adoption rate by race matrix
rgs          = zeros(2,1);                              % Create a vector of racial group sizes

[I J]              = find(nrig(:,2)==1);                % Find whites
rgs(1)             = size(I,1);                         % Count the number of whites
adopt_race(1,:)    = sum(adopt_ind(I,:));               % Count the number of white adopters over time
rate_race(1,:)     = adopt_race(1,:)./rgs(1);           % Compute the rate of adopters over time

[I J]              = find(nrig(:,2)==2);                % Find blacks
rgs(2)             = size(I,1);                         % Count the number of blacks
adopt_race(2,:)    = sum(adopt_ind(I,:));               % Count the number of black adopters over time
rate_race(2,:)     = adopt_race(2,:)./rgs(2);           % Compute the rate of black adopters over time

return;
