program map parameter(mt=9,ndy=14,ny=2009,sc=0.07) c******************************************** c Hinich spatial ideological map program Version 9-14-2009 c Method uses a set of scores which depend on Euclidean distances between c each respondent's ideal ideological point and the candidates positions c in the latent ideological space. c Reference: The Spatial Theory of Voting by Enelow & Hinich c Cambridge University Press (1984) & M. J. Hinich, (2005),"A New Statistical c Multidimensional Unfolding Method," Communication in Statistics - Theory c and Methods, 34: 2299-2310 c c================================================= c c Lines for the control file map.cnl c IMPORTANT: Comments are put after the symbol * c It is vital that the symbol * be used to indicate the start of the comment. c This delimeter must be placed in a column < = 80. Do not exceed the 80 c column field for control statements. The parameters and control character c statements must be to the left of *, the delimeter. c Each line in map.cnl must have the symbol * in each line. c c===================================== c c Use the map-temp.cnl in the Map folder on my webpage as a template for the c control file map.cnl. c Use the keywords followed by an = sign and the setting. Some settings are c character strings while others are real numbers. c c===================================== c c The following parameters must be set in up to 11 lines in map.cnl. They can be c set in any order providing that the keywords are used and the values are c placed after the equals symbol = . c For example: missing data=% where % is wild card symbol for the missing c data numerical flag such as 99. Another example is rotation=7.0 to rotate c the map seven degrees. c c IMPORTANT: Place the symbols # * in each line after these parameters c are set. The program searches for the symbol # to delimite the parameter c reads followed by the read delimier symbol * after #. c The program then reads the comment line 'Candidate or party list' c and then reads the candidate/party names in the list. c c===================================== c c Path & data file name c c Example: \data\map\file name=2004.data */File name c c===================================== c c Reference candidate=%name or reference party=%name c c The program places this candidate/party at the origin of the map or at the c origin scaled from the point set in the cnl file as explained below. c c The column of scores for this candidate (or party) is subtracted from the c scores of the other candidates. Although the theory presented in c Enelow & Hinich (1984) implies that this choice shoould not matter, in c practice in does matter. You should use a candidate/party that you believe c to be in the middle of the voter ideal points. c c===================================== c c Spatial distance model=Euclidean or squared Euclidean c c For the Euclidean distance model the voter i's grade of candidate C is c assumed to be a linear function of the simple Euclidean distance between c the spatial position of candidate C and voter i's ideal point. c In the squared distance model the grade is assumed to be linear in the c square of the simple Euclidean distance. c The program looks for the character string 'squar' on the line. If there c is no 'squar' present then the default choice is Euclidean distance. c c===================================== c c Political space dimension=1 or 2 (or the words one or two) c c The program looks for the number or the word. c c===================================== c c Voter id=yes or no c c Yes will make the program read a respondent lable from column one of the c data file. This column of characters must be in the data file if the c switch is set =yes. c c===================================== c c Candidate valence dimension=yes or no c c Type yes if your data includes a valence dimension score. A valence c dimension score is a measure of a positive candidate/party c characteristic such as integrity or competence to govern. c If yes is set then a values must be placed as shown in map-temp.cnl c after the candidate (or party) name. c c===================================== c c Output files=yes or no c c Yes will open files Choices.graph & Voters.graph which are data files c to be used by any graphics package. c Choices.graph has the candidate (or party) positions estimated by map. c Voters.graph has the estimated voter ideal points. c c===================================== c c Mean set as follows=yes c Yes will set the mean to the scaling of the coordinates (mean-x,mean-y) c where mean-x & mean-y are set in Map.cnl as follows using wild card c symbols for the two real numbers: c c mean-x=%x mean-y=%y c c Example: mean-x=17.0 mean-y=-7.0 */ set mean coordinates c c The mean coordinates are scaled by the estimated standard deviations (references). c c===================================== c c Set rotation angle=yes c Yes will force a rotation whose angle in degrees is set in Map.cnl as c follows where the wild card % is the angle in degrees: c rotation angle=% c c===================================== c c Reflect axis=one (or two or one & two) to reflect the candidate/party map c through the axis. You can type either 'axis' or axes". c Two examples: c reflect axes=one & two */ Reflect axes c reflect axes=1 & 2 */ Reflect axes c Leave the line blank up to the delimitor * for no reflections or use c reflect axes= */ Reflect axes c The default is a blank line - no reflections c c===================================== c c Missing data=% c where % is the missing data flag. Example: missing data=99 c c Any row that has one or more %'s will be not be read. That respondent's c scores will not be used. c c============================================ c c Map method=original or new c c The default should be 'original' until I finish testing the maximum c likelihood algorithm in the new approach described in the paper c "A New Statistical Multidimensional Unfolding Method," Communication c in Statistics 34, 2299-2310, (2005) c c IMPORTANT: Do not type the symbol ' in the string. The symbol ' is used c to set off the character string from the instruction sentences above. c c===================================== c c The following lines are NOT read by the program if 'issue flags columns=0' c is set on the second line of the data file (see below). c If columns=0 is typed on the second line of the data file then the following c should be omitted from the cnl file: c c issue means flags-col1=%,%:%,% issue/type flag groups=%ng1 c issue means flags-col2=%,%:% issue/type flag groups=%ng2 c ... c issue means flags-col19=%,%:%,%:%,% issue/type flag groups=%ng19 c issue means flags-col20=% issue/type flag groups=%ng20 c c where %ng is the number of flag groups for each whose median idea points are c to be calculated where %ng=1 or 2 or 3. The flags are NON-NEGATIVE INTEGERS. c c Use as many columns as you need up to 20 & up to 3 groups of 1 or 2 flags c IMPORTANT: These flags are integers from 0 to 9. c c The integer flags are used to select respondents with a given issue score or c a certain characteristic determined by up to 20 columns of numbers to the c right of the last grade score column. c c The number of flag columns is set in the second line of the data file. c If 'issue flags columns=0' in the second line of the data file then only the c median of the estimated ideal points is computed. c c The flags are used to select to find the median ideal point of respondents c with certain characteristics determined by one or several columns of c numbers to the right of the last grade score column. c c Type up to 3 group of one or two integers from 0 - 9 separated by commas c after the symbol = for each column to select respondents to compute their median c ideal points. Separate the integers by commas in each group and separate the c groups by the symbol : but there should be a blank at the end of the flags. c c Example 1: c issue/type flag groups=3 c issue means flags-Col1=1,2:5:6,7 c c If the kth respondent's plot flag type in the column equals one of the flags c on the group then that respondent's ideal point estimate is used in the c average. For example 1 and the issue/type column 1 you want to find the 3 c median ideal points of respondents whose responses are 1 or 2 , 5, & 6 or 7. c c Example 2: c issue/type flag groups=1 c issue means flags-Col2=7 c c For example 2 and the issue/type column 2 you want to find the mean ideal c point of respondents whose response is 7. c c===================================== c c Respondent map selection flags=%,%,%,%,% (1 <= 10 integer flags 0 - 9) c c The flags are used to select respondents with certain characteristics c determined by one column of numbers to the right of the last issue median flag c column. If there are no issues median flag columns the map flag column is to c the right of the last grade score column. c c This line is NOT read if 'selection flags=0' is set on the second line c of the data file (see below). c c Type up to 10 integers separated by commas after the symbol = to c select the respondents who are used to estimate the political map. c The word 'respondent' can be omitted. c c Example: map selection flags=1,2,7 (no comma after the last flag) c Separate the numbers by commas with no comma after the last integers. c c The selection method is the same as for the input of plot flag columns c but the flags numbers can be different. c To omit a set of respondents whose flag is say 1 do not use 1 for any c of the flags. c c If 'selection flags=0' in the second line of the data file then all the c respondent's scores are used. c If so then this line can be omitted from Map.cnl c c If there is no column of respondent map flags in the data file then you c must type 'selection flags=0' on the second line of the data file in order c to stop the program from trying to find that column. c c===================================== c c # */ Candidate or party list with valence scores if available c Names of the m candidates/parties underneath the above title c c If there are valence scores type them as v=% on the same line as the c candidate/party name as for example v=1 for candidate 1. c The larger the valence score for a candidate, the higher the grade the c candidate gets from each respondent. c c Type 'skip' or 'no' to the right of a candidate/party you do not want c to use. The program will rearrange the data for you. c c Type 'vote' to the right of a candidate if you want to estimate the c proportion of the respondents who will vote for that candidate/party. c Type at least ONE vote in the list. c c IMPORTANT: Place the symbol * after each name. The program searches for c the symbol * to find the end of the name, valence & vote input. c c Also place the symbols # * in a line after these names are read. c The program searches for the symbol # to find the end of the list. c The delimiter # ends the input of the candidate or party list. c c============================================ c c Bootstrapping the MAP generated by the new method c c If the new method is used then you can bootstrap the data to obtain bootstrapped c standard errors. The resampled bootstrap is controlled as follows: c c Bootstraps=% c c The data is resampled % times with replacement to obtain an estimate of the c sampling properties of the estimates of the mean ideal point, rotation, and c standard deviations of the ideal points using the standard bootstrap method. c No bootstrapping is done for the old method and the line for bootstraps can c be omitted from Map.cnl. I put this bootstrap line at the bottom of Map.cnl c when I use the original method to avoid pasting it in when I switch to the c new method. c c============================================ c c The respondent ideal points are estimated after the candidate map is c estimated. The ideal point subroutine uses a three level approach: c c 1. If a respondent gives the highest possible score to one candidate c the ideal point of the respondent is placed at the estimated position c of that candidate. If the respondent gives the highest score to c several candidates, the respondent's ideal point is placed at the c centroid of those candidate's estimated positions. c The number of candidates with the top score is written out to the c right of the ideal point coordinates. c c 2. If no candidate gets the highest score but one gets the lowest c possible score, the ideal point is determined from the low score and c the candidate who gets the highest score. c For example suppose that the highest score possible is 100 and the c lowest score possible is 0. c Suppose that the respondent gives candidate k a 0 and candidate j a c score of 90. Then his ideal point is a prorated point on the line c segment that connects the positions of candidates j & k. c Only one candidate can have the highest score for any respondent in c this proration class. c The integer 11 is written to the right of the ideal point coordinates. c c 3. Otherwise a least squares fit is made to estimate the ideal point. c The least squares fit is the default. c The integer 21 is written to the right of the ideal point coordinates. c c************************************************************************* c c Organization of the data file structure c c========================================================================== c Line 1: Name of the data c Line 2: respondents=%n, candidates=%m, issue flags columns=%k, selection c flags=0/1, format=(%) c c %n - number of respondents (or regions) c %m - number of candidates (or parties c %k - number of issue means flag columns. If 0 then there are no such flag c columns to be read and only the median value of the estimated ideal points of c all the respondents with no missing data values is computed. c %0/1 - 0 or 1 for a map flag column. If 0 then there is no map flag column c c Line 3: Names of the candidates/parties over the grade columns & the flag c columns. The order of the names must match the order of the names set in c map.cnl. The program checks the order on the data file against the order in c the cnl file. c c The data file has n rows after the three line header. c c If there is respondent identification code it should be in Column 1. c The respondent id is read as a character string from this column if c the switch id=yes is set in map.cnl. c If there is no respondent id code column then the first column is c the scores for the first candidate/party. c c The candidate/party grades are in m adjacent columns of the same c width. The grades must be right justified. If is vital that the c the grades columns have the same width. c c A missing grade is coded as % where % is chosen by the user. c Any respondent who has at least one % is NOT USED. c c RESPONDENT SELECTION FLAG COLUMNS: c c The plot flag columns to the right of the score columns if flags are used. c There must be a column symbol over each column as is the case for the c candidates. c c The map flag column must be the LAST column if it exists. c The plot flag columns contains numbers that classify the respondents c whose ideal points whose means are calculated. c c Format - enclosed in ( ) for the scores data array in 50 characters or less. c Suppose that there are ten candidates with valence scores and a c respondent five character id column and two respondent flag columns c whose maximum value is 99. Assume that the two flag columns have a c one space separation. Also assume that the candidate grades have a c maximum value of 100 and a minimum value of 0 and that all columns c are separated by one space. c c The format on line 2 of the data file must be of the form c (a5,20(1x,f3.0),2(1x,i1)). The a5 format calls for a five character c respondent id. The f3.0 calls for a number of at most three spaces c (such as 100) with nothing to the right of the decimal point or no c decimal point. The 1x is a one space separation between the numbers c on each row. c If the grades are between 0 & 10 use (a5,20(1x,f2.0),2(1x,i1)). c Another example is (10(f7.1),i1,i1) if one is using % aggregate data for c ten candidates or parties and two adjacent integer flag columns. c c A format (10(f7.1),i1,i1) for the data structure requires that each c row contain ten real numbers of the form -12345.6 or 12345.6 lined up c with one value to the right of the decimal and an 11th column with two c adjacent integers from 0 to 9, such as 17. It is better to use 1x c one space separators in order to read the data file. c c If the first column contains lables for the respondents that are only c three characters wide then use a3 instead of a5. c c If there are ten candidates and two flag columns and grades from 0 - 100 c then use (a5,10(1x,f3.0),2(1x,i1)). c c If there is only one flag column then use (a5,10(1x,i3),1x,i1) c for ten candidate grades. c c Any book or manual on Fortran F77 or F90 provides examples of f and i formats. c Always use the same format for each column. c===================================== c Number of candidates required to make a spatial map for dim = 2 c Valence scores for candidates => 8 No valence score in the data => 7 c============================================ allocatable::bf,can,cand,check,flag,ng,uc,uv,va real va(:),xc(2) integer flag(:,:),flgd(10),ng(:),uc(:),uv(:),vi,pow,dim,cns(2) character*1000 name character*80 buf,fmt character*50 charead character*30 can(:),cand(:),check(:),csub character*25 bf(:) character*20 par character*2 delim,ref character acm,ch,dav logical exs c******************************************** open(9,file='Map-error.out') c Open control and summary output file open(1,file='map.cnl',err=1,iostat=io1,status='old') open(3,file='Map.out',err=2,iostat=io2) c Start time call cpu_time(gs) c------- ip=0 do k=1,70 c Clear buf do n=1,80 buf(n:n)=' ' enddo c Read line read(1,'(a80)') buf ia=index(buf,'#') if(ia>0) then exit endif ig=index(buf,'*') if(ig==0) then write(3,*)'Place comment marker * after parameters on map.cnl' write(3,'(a80)')buf write(6,*)'Place comment marker * after parameters on map.cnl' write(6,'(a80)')buf write(9,*)'Place comment marker * after parameters on map.cnl' write(9,'(a80)')buf stop endif ig=ig-1 name(ip+1:ip+ig)=buf(:ig) ip=ip+ig enddo c Data file name from control file par='ile name' delim(1:2)='= ' fmt=charead(par,name,delim,ia,ier,3) if(ier>0) then write(3,*)'Error in reading the file name from map.cnl' write(3,*)'Error number for read = ',ier write(6,*)'Error in reading the file name from map.cnl' write(6,*)'Error number for read = ',ier write(9,*)'Error in reading the file name from map.cnl' write(9,*)'Error number for read = ',ier stop endif c------- close(2) c Inquire if file exists inquire(file=fmt,exist=exs) if(.not.exs) then write(3,*)'Input file ',fmt,' does not exist' write(6,*)'Input file ',fmt,' does not exist' write(9,*)'Input file ',fmt,' does not exist' stop endif open(2,file=fmt,err=3,iostat=io3,status='old') rewind(2) do k=1,80 buf(k:k)=' ' enddo c Read & write file identity call readmap(buf,n,m,kis,mpflg,fmt,2,3) c Trap fmt length i1=index(buf,' ') if(i1>=51) then write(3,'(/4x,''ERROR - format has more than 50 characters'')') write(6,'(/4x,''ERROR - format has more than 50 characters'')') write(9,'(/4x,''ERROR - format has more than 50 characters'')') stop endif c------- if(kis==0) then c No issue flags ks=0 elseif(kis>0) then ks=kis-1 else write(3,'(/4x,''Error in read of the number of issue median f *lags''/)') write(6,'(/4x,''Error in read of the number of issue median f *lags''/)') write(9,'(/4x,''Error in read of the number of issue median f *lags''/)') stop endif c------- c Allocate space m1=m-1 m2=m1+kis allocate(bf(0:ks),can(0:m1),cand(0:m1),check(0:m2),uc(0:m1),uv(0:m *1),va(0:m1),flag(0:ks,0:5),ng(0:ks),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate work space for arrays in Map' stop endif if(kis==0) then c Set flag if kis=0 do k=0,5 flag(0,k)=-1 enddo endif write(3,'(a80)')buf c Spatial distance model ier=0 par='odel' c Euclidean distance pow=2 buf=charead(par,name,delim,ia,ier,3) i1=index(buf,'squar') if(i1>0) then pow=1 endif if(ier>0) then write(3,*)'Distance model switch is wrong' write(3,*)'Check spelling in map.cnl' write(6,*)'Distance model switch is wrong' write(6,*)'Check spelling in map.cnl' write(9,*)'Distance model switch is wrong' write(9,*)'Check spelling in map.cnl' stop endif c Spatial dimension par='dimension' buf=charead(par,name,delim,ia,ier,3) i1=index(buf,'one') i1n=index(buf,'1') i2=index(buf,'two') i2n=index(buf,'2') if(i1>0.or.i1n>0) then dim=1 elseif(i2>0.or.i2n>0) then dim=2 else write(3,*)'Political space dimension is not one or two' write(3,*)'Check spelling in map.cnl' write(6,*)'Political space dimension is not one or two' write(6,*)'Check spelling in map.cnl' write(9,*)'Political space dimension is not one or two' write(9,*)'Check spelling in map.cnl' stop endif c Valence dimension, voter id vi=0 par='nce dimension' buf=charead(par,name,delim,ia,ier,3) if(ier==0) then i1=index(buf,'yes') if(i1>0) then vi=1 endif endif if(ier>0) then write(3,*)'Error in read of valence dimension switch ',ier write(6,*)'Error in read of valence dimension switch ',ier write(9,*)'Error in read of valence dimension switch ',ier stop endif c Voter id dav='n' par='oter id' buf=charead(par,name,delim,ia,ier,3) if(ier>0) then write(3,*)'Error in read of voter id ',ier write(6,*)'Error in read of voter id ',ier write(9,*)'Error in read of voter id ',ier stop endif i1=index(buf,'yes') if(i1>0) then dav='y' endif c------- c Check for switches for origin or rotation contraints cns(1)=0 cns(2)=0 r=0.0 par='ordinates' buf=charead(par,name,delim,ia,ier,3) i1=index(buf,'yes') i2=index(buf,'YES') if(i1>0.or.i2>0) then c Constrain origin to (origin-x,origin-y) cns(1)=1 par='ean-x' xc(1)=rnumread(par,name,delim,ia,ier,io) par='ean-y' xc(2)=rnumread(par,name,delim,ia,ier,io) endif par='otation' buf=charead(par,name,delim,ia,ier,3) i1=index(buf,'yes') i2=index(buf,'YES') if(i1>0.or.i2>0) then c Constrain rotation to r cns(2)=2 c Rotation angle par='ngle' r=rnumread(par,name,delim,ia,ier,io) endif c Reference candidate or party par='andidate' csub=charead(par,name,delim,ia,ier,3) if(ier==1) then par='arty' csub=charead(par,name,delim,ia,ier,3) endif c Trap scatter' if(sc<0.) then write(3,*)'Scale ',sc,' < 0' write(6,*)'Scale ',sc,' < 0' stop elseif(sc>.5) then write(3,*)'Scale ',sc,' > 0.5' write(6,*)'Scale ',sc,' > 0.5' stop endif c Iteration difference for max likelihood search del=0.001 c------ c Reflect axis ia=0 ref='nn' par='eflect axes' buf=charead(par,name,delim,ia,ier,3) if(ier==1) then par='eflect axis' buf=charead(par,name,delim,ia,ier,3) endif if(ia>0) then buf(:12)=name(ia:ia+12) i1=index(buf(:12),'one') i2=index(buf(:12),'1') if(i1>0.or.i2>0) then ref(1:1)='y' endif i1=index(buf(:12),'two') i2=index(buf(:12),'2') if(i1>0 .or. i2>0) then ref(2:2)='y' endif endif c------- c Missing data flag par='ing data' smis=rnumread(par,name,delim,ia,ier,io) if(ier>0) then mis=numread(par,name,delim,ia,ier,io) if(ier>0) then write(3,*)'Error in read of missing data flag' write(6,*)'Error in read of missing data flag' stop endif endif smis=float(mis) c Map method par='method' buf=charead(par,name,delim,ia,ier,3) ia=index(buf,'ginal') ib=index(buf,'new') if(ia>0) then ch='c' elseif(ib>0) then ch='n' else write(3,*)'Error in read of map method' write(6,*)'Error in read of map method' stop endif if(ch=='n') then c Resamples par='straps' nstr=numread(par,name,delim,ia,ier,3) if(ier>0) then write(3,*)'Error in read of bootstraps' write(6,*)'Error in read of bootstraps' stop endif if(nstr<0) then write(3,*) 'Error in read of bootstraps - ',nstr write(6,*) 'Error in read of bootstraps - ',nstr stop elseif(nstr==0) then nstr=1 endif endif c Respondent map selection flags if(mpflg==1) then i1=index(name,'ion flags') if(i1==0) then write(3,'(/4x,''The selection flags line is not in Map.cnl''/)') write(6,'(/4x,''The selection flags line is not in Map.cnl''/)') write(9,'(/4x,''The selection flags line is not in Map.cnl''/)') stop endif if(i1>0) then par='ion flags' buf=charead(par,name,delim,ia,ier,3) ng(0)=1 call flagread(1,ng,buf,nfld,flgd,3) endif endif if(mpflg==0) then c Use all respondents with -1 do k=1,10 flgd(k)=-1 enddo endif c Issue means flags for ideal points i1=index(name,'eans flag') c Trap absense of flag lines in cnl if(i1==0.and.kis>0) then write(3,'(/4x,''The issue means flag list is not in the control f *ile''/)') write(6,'(/4x,''The issue means flag list is not in the control f *ile''/)') write(9,'(/4x,''The issue means flag list is not in the control f *ile''/)') stop endif if(kis>0) then do krw=0,ks c No. of issue/type flag groups if(krw<9) then acm=achar(krw+49) par='ol' par(3:3)=acm bf(krw)=charead(par,name,delim,ia,ier,3) if(ier>0) then write(3,*)'Error in read of issue/type group no. row = ',krw+1 write(6,*)'Error in read of issue/type group no. row = ',krw+1 write(9,*)'Error in read of issue/type group no. row = ',krw+1 stop endif ik=index(name,par(:3)) buf(:20)=name(ik+9:ik+29) ik=index(buf(:20),'=') ng(krw)=ichar(buf(ik+1:ik+1))-48 endif if(krw>=9.and.krw<19) then acm=achar(krw+39) par='ol1' par(4:4)=acm bf(krw)=charead(par,name,delim,ia,ier,3) if(ier>0) then write(3,*)'Error in read of issue/type group no. row = ',krw+1 write(6,*)'Error in read of issue/type group no. row = ',krw+1 write(9,*)'Error in read of issue/type group no. row = ',krw+1 stop endif ik=index(name,par(:4)) buf(:20)=name(ik+9:ik+29) ik=index(buf(:20),'=') ng(krw)=ichar(buf(ik+1:ik+1))-48 endif enddo call flagread(kis,ng,bf,flag,3) endif c Write name of program if(ch=='n') then write(3,'(2x,''The Hinich Method for Spatial Analysis of Politica *l Spaces'')') endif if(ch=='c') then write(3,'(2x,''The Cahoon-Hinich Method for Spatial Analysis of P *olitical Spaces'')') endif c------- write(3,'(/7x,'' Number of respondents = '',i7)')n write(3,'(7x,'' Number of candidates in the data file = '',i3/)') *m if(pow==2) then write(3,'(2x,''Scores are Linear in Euclidean Distance'')') else write(3,'(2x,''Scores are Linear in Squared Euclidean Distance' *')') endif if(dim==2) then write(3,'(/7x,''Political Space is Two Dimensional'')') else write(3,'(/7x,''Political Space is One Dimensional'')') endif if(ref(1:1)=='y' .and. ref(2:2)=='n') then write(3,'(/7x,''Map is reflected through axis one'')') endif if(ref(1:1)=='n' .and. ref(2:2)=='y') then write(3,'(/7x,''Map is reflected through axis two'')') endif if(ref(1:1)=='y' .and. ref(2:2)=='y') then write(3,'(/7x,''Map is reflected through axes one & two'')') endif if(cns(1)==1) then write(3,'(/7x,''Unscaled mean ideal point is set = ('',f5.2,'' , *'',f5.2,'')'')')xc endif if(cns(2)==2) then write(3,'(/7x,''Rotation set to angle '',f7.1)')r endif c------- if(dav=='y') then ia=index(fmt,'a') ib=index(fmt,'(a') if(ia==0) then write(3,'('' There is no a in the format on line 2 of the data f *ile to read the voter id lable''/)') write(3,'(a50)')fmt write(3,'(/'' Correct the format for the first column with the v *oter indentification values'')') write(6,'('' There is no a in the format on line 2 of the data f *ile to read the voter id lable''/)') write(6,'(a50)')fmt write(6,'(/'' Correct the format for the first column with the v *oter indentification values'')') write(9,'('' There is no a in the format on line 2 of the data f *ile to read the voter id lable''/)') write(9,'(a50)')fmt write(9,'(/'' Correct the format for the first column with the v *oter indentification values'')') stop elseif(ib==0) then write(3,'('' The a character format in the format on line 2 of * the data file is not of the form (a##,''/)') write(3,'(a50)')fmt write(3,'(/'' Correct the a format'')') write(6,'('' The a character format in the format on line 2 of * the data file is not of the form (a##,''/)') write(6,'(a50)')fmt write(6,'(/'' Correct the a format'')') write(9,'('' The a character format in the format on line 2 of * the data file is not of the form (a##,''/)') write(9,'(a50)')fmt write(9,'(/'' Correct the a format'')') stop endif endif c Output ideal points to file 7 open(4,file='Choices.graph',err=4,iostat=io4) open(7,file='Voters.graph',err=5,iostat=io5) call vot(can,cand,uc,uv,va,flgd,m,n,pow,dim,sc,csub,dav,cns,vi,fmt *,ref,r,smis,del,xc,check,ch,flag,kis,nfld,ng,nstr,3) write(3,*) c+++++++++++++++++++++++++++++++++++++++++++ c Time and date of run & run time call head(mt,ndy,ny,gs,rt,3) c+++++++++++++++++++++++++++++++++++++++++++ deallocate(bf,can,cand,check,flag,ng,uc,uv,va) stop c------- 1 write(6,*)'**** Error on opening control file ',io1 stop 2 write(6,*)'**** Error on opening output file ',io2 stop 3 write(3,*)'**** Error on opening data file ',io3 stop 4 write(3,*)'**** Error on opening candidates/parties file ',io4 stop 5 write(3,*)'**** Error on opening voter ideal point file ',io5 stop end c c=========================================== c subroutine vot(can,cand,uc,uv,va,flgd,m,n,pow,dim,sc,csub,dav,cns, *vi,fmt,ref,r,smis,del,xc,check,ch,flag,kis,nfld,ng,nstr,io) c******************************************** c m - no. candidates in list, mu - no. used, nstr - resamples, del - iteration mesh c nfld - flags in flgd, kis - issue median columns c cns=1 - constrain mean, cns=2, r - rotation c============================================ allocatable::a,cr,cov,d,dbar,e,er,f,s,v,vs,vd,w,x real x(:),cr(:),va(m),vs(:),xc(2) real*8 a(:),cov(:),d(:),dbar(:),v(:),er(:),f(:),e(:),s(:),w(:) integer uc(m),uv(m),flag(kis,6),flgd(10),ng(kis),vi,p *ow,dim,psub,cns(2) character*700 name character*80 buf character*50 charead,fmt,fmo character*30 csub character*30 can(m),cand(m),check(m+kis) character*20 par character*5 vd(:) character*2 delim,ref character ch,dav intent(in) flgd,m,n,pow,dim,csub,cns,vi,fmt,ref,sc,r,smis,del,xc,c *h,flag,kis,nfld,ng,nstr,io c******************************************** c Read candidate names mu=0 psub=0 nsum=0 vmax=0.0 do k=1,m read(1,'(80a)') name c Check for blanks ia=index(name(1:1),' ') ib=index(name(2:2),' ') ic=index(name(3:3),' ') if(ia>0 .and. ib>0 .and. ic>0) then write(io,'('' There are 3 blank columns in the candidate/party n *ame in Map.cnl'')') write(io,'(a77)')name(:77) write(6,'('' There are 3 blank columns in the candidate/party na *me in Map.cnl'')') write(6,'(a77)')name(:77) write(9,'('' There are 3 blank columns in the candidate/party na *me in Map.cnl'')') write(9,'(a77)')name(:77) stop endif ia=index(name,'#') if(ia>0) then if(k-1/=m) then write(io,'(/'' ERROR''/)') write(io,'('' The number of candidates '',i2,'' set in the data * file is not the same as the'',/,'' number of candidates in map.cn *l'')')m write(6,'(/'' ERROR''/)') write(6,'('' The number of candidates '',i2,'' set in the data *file is not the same as the'',/,'' number of candidates in map.cn *l'')')m write(9,'(/'' ERROR''/)') write(9,'('' The number of candidates '',i2,'' set in the data *file is not the same as the'',/,'' number of candidates in map.cn *l'')')m stop endif exit endif ig=index(name,'*') if(ig==0) then write(io,*)'Place comment marker * after candidate/party names a *nd vote switches in the list' write(6,*)'Place comment marker * after candidate/party names an *d vote switches in the list' write(9,*)'Place comment marker * after candidate/party names an *d vote switches in the list' write(io,'(/1x,a77)')name(:77) write(6,'(/1x,a77)')name(:77) write(9,'(/1x,a77)')name(:77) stop endif va(k)=0.0 uc(k)=0 ia=index(name,'=') if(ia>0) then cand(k)=name(:ia-3) else c Find blank after character name ib=index(name,' ') cand(k)=name(:ib) endif i1=index(name(4:),' no ') i2=index(name(4:),' skip ') if((i1>0.or.i2>0).and.cand(k)==csub) then write(io,'(2x,''ERROR - you are skipping the reference candidate * '',i2,1x,a30)')k,csub write(6,'(2x,''ERROR - you are skipping the reference candidate *'',i2,1x,a30)')k,csub write(9,'(2x,''ERROR - you are skipping the reference candidate *'',i2,1x,a30)')k,csub stop endif if(i1==0.and.i2==0) then c Candidate/party k is mapped uc(k)=1 mu=mu+1 can(mu)=cand(k) if(can(mu)==csub) then psub=mu endif c Find vote count flag ia=index(name,'vote') uv(mu)=0 if(ia>0) then c Compute candidate vote count uv(mu)=1 nsum=nsum+1 endif delim='=' if(vi==1) then c Find valence value if v= is present par='v' ve=rnumread(par,name,delim,ia,ier,io) va(mu)=ve if(ier>0) then write(6,*)'Error in read of valence ',k write(io,*)'Error in read of valence ',k stop endif if(ve>vmax) then vmax=ve endif endif endif enddo c End of candidate loop do k=1,mu c Highest valence is 0 va(k)=vmax-va(k) enddo c------- c Trap psub if(psub==0) then write(io,*)'Check spelling of the name of the reference candidate *' write(io,*) csub write(6,*)'Check spelling of the name of the reference candidate' write(6,*)csub stop endif c Trap nsum if(nsum==0) then write(io,'(/'' No candidate or party vote is being calculated'')' *) write(io,'('' Type vote to the right of a candidate or party to g *et a test vote'')') write(6,'(/'' No candidate or party vote is being calculated'')') write(6,'('' Type vote to the right of a candidate or party to ge *t a test vote'')') stop endif c Read candidate name line read(2,'(a700)') name call lineread(name,check,m+kis,io) do kc=1,m if(cand(kc)/=check(kc)) then write(io,'(/'' ERROR''/)') write(io,'('' Candidate-Party Name '',a30,/'' is not the same as * the name over column '',a30)') cand(kc),check(kc) write(io,'(/7x,''Verify that the number of candidates '',i2,'' s *et in the data file'',/,'' is the same as the number of candidates * in map.cnl''/)')m write(io,'(2x,''Verify that there is at least one space between * the candidate names in the 3rd line of the data file'')') write(6,'(/'' ERROR''/)') write(6,'('' Candidate-Party Name '',a30,/'' is not the same as *the name over column '',a30)') cand(kc),check(kc) write(6,'(/7x,''Verify that the number of candidates '',i2,'' se *t in the data file'',/,'' is the same as the number of candidates *in map.cnl''/)')m write(6,'(2x,''Verify that there is at least one space between t *he candidate names in the 3rd line of the data file'')') stop endif enddo c------- if(vi==1) then write(io,'(/7x,''Candidate/Party Valence Scores Used'')') endif write(io,'(/7x,a20,'' scores are subtracted from the others'')') c *an(psub) if(cns(1)==0) then write(io,'(/7x,''The position of this candidate or party is set a *t the origin'')') endif if(cns(1)==1) then write(io,'(/7x,''The position of this candidate or party is set a *t the mean'')') endif write(io,'(/2x,67(''=''))') c------- do k=1,50 fmo(k:k)=' ' enddo if(mu==m .and. dav=='n') then c No respondent id in first column & all candidates are mapped ia=index(fmt,'a') if(ia>0) then c Replace a# by #x, to remove id character read ib=index(fmt(ia+1:),',') fmo(1:1)='(' fmo(2:ib)=fmt(ia+1:ia+ib-1) nt=ib+1 fmo(nt:nt)='x' fmo(nt+1:)=fmt(ia+ib:) else c No a in fmt, shift to fmo fmo(:50)=fmt(:50) endif endif c------- c Respondent id in first column & all candidates are mapped if(mu==m .and. dav=='y') then c Shift to fmo with voter id format fmo(:50)=fmt(:50) endif c No respondent id in first column & some candidates are skipped if(mu0) then c Replace a# by #x, to remove id character read ib=index(fmt(ia+1:),',') fmo(1:1)='(' fmo(2:ib)=fmt(ia+1:ia+ib-1) nt=ib+1 fmo(nt:nt)='x' fmo(nt+1:)=fmt(ia+ib:) else c No a in fmt, shift to fmo fmo(:50)=fmt(:50) endif endif c Respondent id in first column & some candidates are skipped if(mu=0) then write(io,'(/2x,i5,'' respondents for estimating the map were sele *cted using the flags:'',:,10(i2,1x)/)')nu,(flgd(k),k=1,nfld) else write(io,'(/7x,i7,'' respondents were used to map the space'')' *)nu endif c------- ns=nu-1 ka=kis*12-1 allocate(c(0:ns,0:mn),h(0:ns,0:ka),us(0:mn),yc(0:3*mn),stat=ibad) if(ibad/=0) then write(io,*)'Unable to allocate space for more arrays in top' stop endif c Rearrange data column x =>c & valence va =>vs call reo(x,c,dbar,can,cand,iw,va,vs,us,uv,e,a,m,mu,n,nu,pow,psub,s *max,smin,smis,flgd,vi,vmax,kw,io) c------- c Covariance matrix of adjusted scores call covm(c,cov,dbar,nu,mn,io) if(ch=='n') then call new(can,cand,a,c,cov,cr,d,dbar,e,er,f,h,s,uc,uv,v,va,vs,vd,w *,x,flgd,m,mu,mn,n,nu,pow,dim,iw,psub,cns,dav,vi,fmo,ref,sc,r,smis, *nstr,del,smax,smin,vmax,xc,med,check,rflg,flag,kis,nfld,ng,nfg,io) endif if(ch=='c') then call geter(cov,v,er,cosn,a,d,f,e,w,mu,mn,nu,dim,io) c Factor loadings =>v if(pow==2) then c Error variances for Euclidean distance do k=1,mn db=dble(dbar(k)+dbar(mu)) cs=db*db-er(k)/2.d0 c Trap negative cs if(cs<=0.d0) then er(k)=er(k)/(4.d0*db) else er(k)=db-dsqrt(cs) endif enddo cs=dble(dbar(mu)*dbar(mu))-er(mu)/2.d0 if(cs<=0.d0) then er(mu)=er(mu)/(4.d0*dble(dbar(mu))) else er(mu)=dble(dbar(mu))-dsqrt(cs) endif c Adjust means for pow=2 do k=1,mn e(k)=dbar(k)+sngl(er(mu)-er(k)) enddo endif if(pow==1) then do k=1,mn e(k)=dbar(k) enddo endif c Error standard deviations write(io,'(/7x,''Specific Error Standard Deviations''/)') write(io,'(2x,a20,2x,f15.4)') (cand(k),dsqrt(er(k)),k=1,mu) call rotate(v,yc,cosn,dim,mn,io) write(io,'(/7x,''OLS Estimates of Remaining Parameters''/)') if(dim==2) then call v2(mu,yc,cr,v,d,e,w,er,vi,va,wv,cns,r,r2,b0,b1,b2,xc,xm,x0, *io) endif if(dim==1) then call v1(mu,yc,cr,v,d,e,w,er,vi,va,wv,cns,b1,r2,io) endif c------- c Write parameters and candidate positions if(dim==2) then write(io,'(/2x,''Ideal Point Standard Deviations = '',2(f14.2,1x *))')b1,b2 if(vi==1) then write(io,'(2x,''Candidate Characteristic Value Weight = '',f14. *2)') wv endif write(io,'(2x,''Rotation = '',f7.1/)')b0 write(io,'(17x,''Positions of the Candidates''/)') do k=1,mu if(ref(1:1)=='y') then c Reflect around axis 1 cr(k,2)=-1.0*cr(k,2) xm(2)=-1.0*xm(2) endif if(ref(2:2)=='y') then c Reflect around axis 2 cr(k,1)=-1.0*cr(k,1) xm(1)=-1.0*xm(1) endif write(io,'(4x,a20,4x,f12.3,2x,f12.3)')cand(k),cr(k,1),cr(k,2) c Output candidates/parties to file 4 write(4,'(a20,1x,2(f12.3,1x))')cand(k),cr(k,1),cr(k,2) enddo write(io,'(8x,''Mean Ideal Point'',4x,f12.3,2x,f12.3)')x0(1),x0( *2) endif if(dim==1) then c Write parameters and candidate positions for one dimension write(io,'(/2x,''Ideal Point Standard Deviation = '',f14.2)')b1 if(vi==1) then write(io,'(2x,''Candidate Characteristic Value Weight = '',f14. *2)')wv endif write(io,'(/17x,''Positions of the Candidates''/)') do k=1,mu write(io,'(4x,a20,4x,f12.3)')cand(k),cr(k,1) enddo endif c-------- c Plot ideal points if(dim==2) then call ideal2(n,nu,mu,c,uv,cr,dbar,va,wv,a,e,iw,y,vd,dav,pow,ref,s *c,io) write(io,'(''Mean of the Estimated Ideal Points'',2x,2(f12.3,2x) */)')y(1),y(2) write(4,'(''Mean Ideal Point'',5x,2(f12.3,1x))')x0(1),x0(2) write(4,'(''xbar'',17x,2(f12.3,1x))')y(1),y(2) call flagmed(c,rflg,m,kis,n,nu,flag,ng,check,h,med,io,4) write(io,'(2x,70(''=''))') c Voting % estimates write(io,'(7x,''% of Ideal Points that are Closest to Each Candi *date or Party''/)') do k=1,mu if(uv(k)==1) then c Vote % write(io,'(a20,1x,f5.1)')can(k),100.*sngl(a(k))/float(nu) endif enddo endif c------- if(dim==1) then call ideal1(nu,mu,c,cr,iw,sc,io) do k=1,nu ncode=nint(x(k,2)) write(7,'(f12.3,1x,i3)')x(k,1),ncode enddo do k=1,mu write(4,'(a20,1x,f12.3)')cand(k),cr(k,1) enddo endif endif c------- deallocate(iw,rflg,h,med,us,yc) return end c c=========================================== c subroutine dat(x,iw,rflg,uc,kis,m,mu,n,nu,smax,smin,fmo,flgd,dav,s *mis,vd,io) c******************************************** c Read data & find max & min of data . c nu - number of respondents used, kis - issue flag columns c flgd - flags for respondents used, smis - missing data flag c Data from all respondents if flgd(1)=-1 c mu - candidates used, m - all candidates, iw - missing data count c uc - candidate use flag, vd - character respondent identification c============================================ allocatable::w,wflg real x(n,mu),w(:),smax,smin integer iw(mu),rflg(n,kis),uc(m),wflg(:),flgd(10),use,kchk character*50 fmo character*5 vd(n),rvd character dav intent(in) kis,m,mu,n,fmo,flgd,dav,smis,uc,io intent(out) x,rflg,iw,nu,smax,smin,vd c******************************************** mn=m-1 allocate(w(0:mn),wflg(0:kis-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate work space for w in dat' stop endif c Output scores for first respondent to check against data file write(io,'(/2x,67(''*''))') write(io,'(2x,''Data for First Respondent to Check Data Read'')') c------- nmp=0 if(flgd(1)/=-1) then c Map flag column nmp=1 endif c Initialize do k=1,mu iw(k)=0 enddo nu=0 noms=0 mis1=0 mis2=0 mis3=0 mis4=0 mis5=0 smax=0. smin=100. c------- c Voter loop do i=1,n nti=-1 c------- if(nmp==0.and.kis==0) then c Input x(i,j) for all i with no flags if(dav=='y') then read(2,fmo) rvd,(w(j),j=0,mn) if(i==1) then write(io,*)rvd,(w(j),j=0,mn) endif endif if(dav/='y') then read(2,fmo) (w(j),j=0,mn) if(i==1) then write(io,*)(w(j),j=0,mn) endif endif endif c------- if(nmp==0.and.kis>0) then c Input x(i,j) for all i =>c(n,j) & issue flags if(dav=='y') then c Read respondent id read(2,fmo) rvd,(w(j),j=0,mn),(wflg(j),j=0,kis-1) if(i==1) then write(io,*)rvd,(w(j),j=0,mn),(wflg(j),j=0,kis-1) endif else c Input x(i,j) for all i with no id read(2,fmo) (w(j),j=0,mn),(wflg(j),j=0,kis-1) if(i==1) then write(io,*)(w(j),j=0,mn),(wflg(j),j=0,kis-1) endif endif endif c------- if(nmp==1.and.kis==0) then c Input data for all i & map flag column with no issue flags if(dav=='y') then read(2,fmo) rvd,(w(j),j=0,mn),nti if(i==1) then write(io,*)rvd,(w(j),j=0,mn),nti endif endif if(dav/='y') then read(2,fmo) (w(j),j=0,mn),nti if(i==1) then write(io,*)(w(j),j=0,mn),nti endif endif endif c------- if(nmp==1.and.kis>0) then c Input data for all i =>c(n,j), issue flags & map flag if(dav=='y') then c Read respondent id read(2,fmo) rvd,(w(j),j=0,mn),(wflg(j),j=0,kis-1),nti if(i==1) then write(io,*)rvd,(w(j),j=0,mn),(wflg(j),j=0,kis-1),nti endif else c Input data for all i with no id read(2,fmo) (w(j),j=0,mn),(wflg(j),j=0,kis-1),nti if(i==1) then write(io,*)(w(j),j=0,mn),(wflg(j),j=0,kis-1),nti endif endif endif c------- knt=0 mis=0 ncmis=-1 do k=1,m if(uc(k)==1) then knt=knt+1 sv=w(k-1) c Check for missing data flag if(sv==smis) then c Missing data flag mis=mis+1 iw(knt)=iw(knt)+1 endif endif enddo if(mis==0) then noms=noms+1 endif if(mis==1) then mis1=mis1+1 endif if(mis==2) then mis2=mis2+1 endif if(mis==3) then mis3=mis3+1 endif if(mis==4) then mis4=mis4+1 endif if(mis==5) then mis5=mis5+1 endif c------- do k=1,m sv=w(k-1) if(mis==0.and.uc(k)==1) then c No missing value & use column if(svsmax) then smax=sv endif endif c End of max&min row sweep enddo c------- c Count selected voters using the flgd flages nuse=0 do ktst=1,10 if(nti==flgd(ktst).or.nti==-1) then nuse=1 endif enddo c Check if scores are different for each other kchk=0 do kw=1,mn if(w(kw)==w(0)) then kchk=kchk+1 endif enddo if(kchk==mn) then nuse=0 endif c------- if(nuse==1) then if(mis==0) then nu=nu+1 use=0 c------- c Move w =>x if only 1 missing & scores are different for each other do kw=1,m if(uc(kw)==1) then c Use column use=use+1 x(nu,use)=w(kw-1) endif enddo c Trap nuse if(use/=mu) then write(io,*)'Number used ',use, '/= candidate used ',mu write(6,*)'Number used ',use, '/= candidate used ',mu stop endif if(kis>0) then c Move wflg =>rflg do kf=1,kis rflg(nu,kf)=wflg(kf-1) enddo endif c------- if(dav=='y') then c Respondent identification =>vd vd(nu)=rvd endif endif endif c End of voter loop enddo c Trap range if(smax<=0) then write(io,*) write(io,*)'Error' write(io,*)'max score ',smax,' not positive' write(6,*) write(6,*)'Error' write(6,*)'max score ',smax,' not positive' stop elseif(smin>=smax) then write(io,*) write(io,*)'Error' write(io,*)'smin ',smin,' > smax ',smax write(6,*) write(6,*)'Error' write(6,*)'smin ',smin,' > smax ',smax stop endif c Trap nu if(nu==0) then write(io,'(/7x,''Error'')') write(io,'(/7x,''Error'')') write(io,*)'Respondents selected by the first flag = 0' write(6,'(/7x,''Error'')') write(6,'(/7x,''Error'')') write(6,*)'Respondents selected by the first flag = 0' stop endif write(io,'(/7x,i6,'' Respondents with No Missing Value'')')noms write(io,'(/7x,i6,'' Respondents with One Missing Value'')')mis1 write(io,'(/7x,i6,'' Respondents with Two Missing Values'')')mis2 write(io,'(/7x,i6,'' Respondents with Three Missing Values'')')mis *3 write(io,'(/7x,i6,'' Respondents with Four Missing Values'')')mis4 write(io,'(/7x,i6,'' Respondents with Five Missing Values'')')mis5 c------- deallocate(w) return end c c=========================================== c subroutine reo(x,c,dbar,can,cand,iw,va,vs,us,uv,e,a,m,mu,n,nu,pow, *psub,smax,smin,smis,flgd,vi,vmax,kw,io) c******************************************** c Rearrange columns x => c . c nu - number of respondents used, ncol - flag columns (0 or 1) c flgd - flags for respondents used, smis - missing data flag c Data from all respondents if flgd(1)=-1 c mu - candidates used, m - all candidates c uv - vote count flag c============================================ integer us(mu),uv(mu),iw(mu),kw(mu),psub,pow,flgd(7),use,vi real x(n,mu),c(nu,mu),va(m),vs(m),smax,smin real*8 e(mu),a(2*mu),dbar(mu),rn character*30 can(m),cand(m),buf intent(in) x,can,iw,m,mu,n,nu,pow,psub,flgd,smax,smin,smis,uv,va,v *i,vmax,io intent(out) c,dbar,us,vs c******************************************** smd=(smax+1-smin)/2 c Voter loop do i=1,nu do j=1,mu sv=x(i,j) if(sv==smis) then c Replace missing value with mid value sv=smd endif c Move x =>c c(i,j)=sv enddo enddo c End of voter loop write(io,'(2x,67(''=''))') c------- c Voter loop for columns shifts if(psub/=mu) then c Rearrange columns of x =>c, candidates =>cand, valence =>vs, uv =>us nb=psub-1 nc=psub+1 buf=can(psub) uh=uv(psub) vh=va(psub) ih=iw(psub) do i=1,nu cst=c(i,psub) c Rearrange nc-mu columns do j=nc,mu c(i,j-1)=c(i,j) if(i==nu) then cand(j-1)=can(j) vs(j-1)=va(j) us(j-1)=uv(j) kw(j-1)=iw(j) if(psub>1) then c Shift can => cand & va =>vs do kc=1,nb cand(kc)=can(kc) us(kc)=uv(kc) vs(kc)=va(kc) kw(kc)=iw(kc) enddo endif endif enddo c(i,mu)=cst enddo c End of column shift loop cand(mu)=buf us(mu)=uh vs(mu)=vh kw(mu)=ih endif if(psub==mu) then do j=1,mu cand(j)=can(j) vs(j)=va(j) kw(j)=iw(j) enddo endif c------- c Second voter loop m2=2*mu do k=1,mu e(k)=0. a(k)=0.d0 a(mu+k)=0.d0 enddo do i=1,nu c Accumulate first two moments of c a(mu)=a(mu)+dble(c(i,mu)) a(m2)=a(m2)+dble(c(i,mu)*c(i,mu)) do j=1,mu-1 c Check smax if(c(i,j)>smax) then write(io,'(/''ERROR''/)') write(io,*)'x(i,j) ',c(i,j),' > max score ',smax,' i & j ',i,j write(io,*) 'This observation exceeds the maximum value' stop endif a(j)=a(j)+dble(c(i,j)) a(mu+j)=a(mu+j)+dble(c(i,j)*c(i,j)) enddo c------- c Transform scores so that 0 is the best score for a candidate if(pow==1) then c(i,mu)=smax-c(i,mu) e(mu)=e(mu)+dble(c(i,mu)) endif c Square if pow=2 & subtract last column if(pow==2) then c(i,mu)=smax-c(i,mu) c(i,mu)=c(i,mu)*c(i,mu) e(mu)=e(mu)+dble(c(i,mu)) endif c Subtract transformed scores from smax =>0 do j=1,mu-1 if(pow==1) then c(i,j)=smax-c(i,j) e(j)=e(j)+dble(c(i,j)) endif if(pow==2) then c(i,j)=smax-c(i,j) c(i,j)=c(i,j)*c(i,j) e(j)=e(j)+dble(c(i,j)) endif c Subtract last squared column from each column c(i,j)=c(i,j)-c(i,mu) enddo c End of second voter loop enddo c------- c Means of smax-score rn=dfloat(nu) dbar(mu)=e(mu)/rn do k=1,mu-1 dbar(k)=sngl(e(k)/rn)-dbar(mu) enddo if(vi==0) then write(io,'(/32x,''Means Standard Deviations Missing Scores'' */)') elseif(vi==1) then write(io,'(/32x,''Means Standard Deviations Valence Scores M *issing Scores''/)') endif do j=1,mu a(j)=a(j)/rn a(mu+j)=a(mu+j)/float(nu-1)-(rn/float(nu-1))*a(j)*a(j) if(vi==0) then write(io,'(1x,a20,2(3x,f12.2),14x,i5)')cand(j),a(j),sqrt(a(mu+j) *),kw(j) elseif(vi==1) then write(io,'(1x,a20,3(3x,f12.2),14x,i5)')cand(j),a(j),sqrt(a(mu+ *j)),vmax-vs(j),kw(j) endif enddo c------- return end c c=========================================== c subroutine covm(x,cov,dbar,nu,mn,io) c******************************************** c Covariance matrix of adjusted scores used c============================================ real x(nu,mn) real*8 cov(mn,mn),dbar(mn),ck1,ck2 intent(in) x,dbar,mn,nu intent(out) cov c******************************************** c Initialize do k1=1,mn do k2=1,mn cov(k1,k2)=0.d0 enddo enddo c Covariance matrix do i=1,nu do k1=1,mn ck1=dble(x(i,k1))-dbar(k1) do k2=k1,mn ck2=dble(x(i,k2))-dbar(k2) cov(k1,k2)=cov(k1,k2)+ck1*ck2 enddo enddo enddo do k1=1,mn do k2=k1,mn cov(k1,k2)=cov(k1,k2)/dfloat(nu) cov(k2,k1)=cov(k1,k2) enddo enddo c------- return end c c=========================================== c subroutine geter(cov,v,er,cosn,a,d,f,e,w,m,mn,nr,dim,io) c************************************************* c Maximum likelihood factor analysis & estimate the error variance of c(m) c m - number of candidates used, nr - number of respondents, mn=m-1 c cov - covariances of adjusted scores, dim - no. of ideological dimensions c v - n x (dim+1) matrix of factor loadings, cosn(1)=cos1 & cosn(2)=sin2 c er - Specific error variances, z1 - min of -ln(L) c================================================= implicit real*8(a-h,o-z) allocatable::c real*8 a(m*m),c(:,:),cov(mn,mn),d(mn*mn),e(m),er(m),f(mn),v(mn,mn) *,w(mn*mn) real cosn(2),r2 integer dim,sw intent(in) cov,m,mn,nr,dim,io intent(out) v,er,cosn c************************************************* allocate(c(0:mn-1,0:mn-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate c in maxle' stop endif nvec=dim+1 c Factor analysis call cfac2(cov,v,nr,mn,nvec,a,d,er,e,f,c,z1,1,ier,io) if(ier==1) then write(io,*) 'Cov is not positive definite in cfac in geter' write(6,*) 'Cov is not positive definite in cfac in geter' stop endif c Dependent variable for first regression do k=1,mn f(k)=1.d0 enddo c First least squares fit do i=1,nvec do j=1,mn v(j,i)=v(j,i)/2.d0 a((i-1)*mn+j)=v(j,i) enddo enddo c------- call regres(f,a,w,e,mn,nvec,0,r2,1,ier,io) c Solve for error variance & rotation paramaters c2=0.d0 do k=1,nvec c2=c2+w(k)*w(k) enddo sc2=dsqrt(c2) er(m)=4./sngl(c2) do k=1,nvec w(k)=w(k)/sc2 enddo c------- if(dim==2) then sin2=sngl(w(1)) cos2=sqrt(1.-sin2*sin2) if(w(2)s from maxle sw=1 call maxle(cov,nu,m,mn,dim,a,d,e,er,f,v,w,z1,del,sw,ier,io) if(pow==2) then c Error variance for Euclidean distance do k=1,mn db=dble(dbar(k)+dbar(m)) cs=db*db-er(k)/2.d0 c Trap negative cs if(cs<=0.d0) then er(k)=er(k)/(4.d0*db) else er(k)=db-dsqrt(cs) endif enddo cs=dble(dbar(m)*dbar(m))-er(m)/2.d0 if(cs<=0.d0) then er(m)=er(m)/(4.d0*dble(dbar(m))) else er(m)=dble(dbar(m))-dsqrt(cs) endif c Adjust means for pow=2 do k=1,mn tm(k-1)=dbar(k)+sngl(er(m)-er(k)) enddo endif if(pow==1) then do k=1,mn tm(k-1)=dbar(k) enddo endif c Set variables for OLS fit & solve for parameters do k=1,mn xr(k-1,0)=v(k,1)*v(k,1) xr(k-1,1)=v(k,2)*v(k,2) xr(k-1,2)=v(k,1)*v(k,2) xr(k-1,3)=v(k,1) xr(k-1,4)=v(k,2) enddo ir=5 c Check for valence weight if(vi>0) then ir=6 do k=1,mn xr(k-1,5)=va(k) enddo endif c OLS fit tm=xr*w int=0 if(kb==1) then c Write r*2 sw=1 else sw=0 endif call regres(tm,xr,w,e,mn,ir,int,r2,sw,ier,io) c Solve for parameters from w call param(w,cns,b1,b2,cosr,sinr,r,ro,xm,io) if(vi>0) then valw=w(6) else valw=0. endif c------- deallocate(tm,xr) return end c c=========================================== c subroutine maxle(cov,n,m,mn,dim,a,d,e,er,f,v,w,z1,del,sw,ier,io) c************************************************* c mn - Number of targets-reference, sw=1 output eigenvalues c dim - Common factors, cov = vv'+sg0*11'+er, sw=2 output all writes c v - matrix of factor loadings, er - Specific variances, z1 - min of -ln(L) c fmin - min of likelihood objective function, sgs - sim(m) estimate c 2nd derivative matrix of S=>d, Lawley's approximation =>a c================================================= implicit real*8(a-h,o-z) allocatable::c real*8 a(m*m),c(:,:),cov(mn,mn),d(mn,mn),v(mn,mn),w(mn,mn),e(m),er *(m),f(mn) real del,sge integer dim,sw,ww,t intent(in) cov,m,mn,n,del,sw,io intent(out) a,d,er,e,ier,v,z1 c************************************************* allocate(c(0:mn-1,0:mn-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate c in maxle' stop endif c Find limit er for loop call cfac2(cov,v,n,mn,dim,a,d,er,e,f,c,z1,0,ier,io) if(ier==1) then write(io,*) 'Cov is not positive definite in first call to cfac' write(6,*) 'Cov is not positive definite in first call to cfac' stop endif sgs=1.d0 do k=1,mn ser=dsqrt(er(k)) if(ser>sgs) then sgs=ser endif enddo ntry=nint(sngl(sgs)/del) ier=0 c Find minimum fac for mle of sg0 t=0 fmin=1.d7 sgs=1.d0 dowhile(t w do k1=1,mn do k2=1,k1 cd=cov(k1,k2)-str*str w(k1,k2)=cd w(k2,k1)=cd enddo enddo c Factor w - no output if sw=0 call cfac2(w,v,n,mn,dim,a,d,er,e,f,c,z1,0,ier,io) if(ier==0) then if(z1d, Lawley's approximation =>a, sw=2 - all writes c================================================= implicit real*8(a-h,o-z) allocatable::indx,w real*8 a(n,n),cov(n,n),d(n,n),s(n,n),v(n,n),er(n),f(n),e(n),w(:),h real dd integer indx(:),sw intent(in) cov,n,nr,nvec,sw,io intent(out) a,d,er,e,v,z1,ier c************************************************* allocate(indx(0:n-1),w(0:n-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate indx in cfac2' stop endif c Max for loop mxt=100 c Initialize z1 & con2 z1=1.0d10 con2=cov(1,1)*7.d-2 kp1=nvec+1 ier=0 c------- c Cov => a & d do i=1,n f(i)=0. er(i)=0.d0 do j=i,n d(i,j)=cov(i,j) enddo enddo inv=0 c w=d**-1*f=0 call dls(d,w,f,e,n,inv,ier,io) c Test if cov is positive definite if(ier==1.and.sw==2) then write(io,*) 'Cov is not positive definite in cfac2' return endif c Initialize er do i=1,n c er=(1-nvec)/(2*n*diag(cov)) er(i)=(1.d0-dfloat(nvec)/(2.d0*dfloat(n)))/d(i,i) if(er(i)<1.d-12) then ier=-1 exit endif enddo if(ier==-1) then return endif c Eigenvalues of (er**-1/2)*cov*(er**-1/2) => e, eigenvectors => v call phi(cov,a,v,e,er,w,n,nvec,io) c 1st derivative of objective S => f call fd(cov,v,e,er,n,nvec,f) c w=a**-1*f call dls(a,w,f,e,n,0,ier,3) if(ier==1.and.sw==2) then write(io,*) 'a is not positive definite' ier=2 return endif c er-w => er do i=1,n er(i)=er(i)-w(i-1) if(er(i)<1.d-10) then er(i)=1.d-10 endif enddo c------- ni=0 z2=0. c Start of MLE loop dowhile (z2<0.99999 .and. ni<=mxt) ni=ni+1 do i=1,n if(er(i)<1.d-10) then er(i)=1.d-10 endif enddo c Eigenvalues of (er**-1/2)*cov*(er**-1/2) => e, eigenvectors => v call phi(cov,a,v,e,er,w,n,nvec,io) if (z1<=0.d0) then z1=1.d-2 endif c 1st derivatives of the log likelihood S => f call fd(cov,v,e,er,n,nvec,f) c 2nd derivative matrix of the log likelihood S => d call sd(v,e,er,n,nvec,f,a,d) c Objective function Sum(e-log(e)-1) z=0.d0 do k=kp1,n z=z+e(k)-dlog(e(k))-1.d0 enddo z2=z1 z1=z c f=>e & d=>s do k=1,n e(k)=f(k) s(k,k)=d(k,k) do j=k+1,n h=d(j,k) s(j,k)=h s(k,j)=h enddo enddo c w=d**-1*f call dludcmp(s,n,n,indx,dd,ier,io) call dlubksb(s,n,n,indx,w) nchk=0 do j=1,n if(dabs(w(j-1))>con2) then nchk=nchk+1 endif if(nchk==1) then do k=1,n w(k-1)=e(k) enddo endif enddo do k=1,n er(k)=er(k)-w(k-1) enddo z2=z1/z2 enddo c End of MLE loop c------- if(ni==mxt) then write(io,'(''Loop terminated at max '',i4,'' iterations'')') mxt ier=3 return endif c Trap negative er's do k=1,n if(er(k)<1.d-10) then er(k)=1.d-10 endif enddo c Lawley's approximation => a call phi(cov,a,v,e,er,w,n,nvec,io) c v=>s do j=1,nvec do i=1,n s(i,j)=v(i,j) enddo enddo c Objective function Sum(e-log(e)-1) z1=0.d0 do k=kp1,n z1=z1+e(k)-dlog(e(k))-1.d0 enddo c sqrt(er(e-1))*v do i=1,nvec do j=1,n ch=e(i)-1.d0 if(ch<=0.d-10) then ch=1.d-10 endif v(j,i)=s(j,i)*dsqrt(er(j)*ch) enddo enddo if(sw==1) then c Output eigenvalues call wrtfac(e,z1,nr,n,nvec,ni,io) endif c------- deallocate(indx,w) return end c c=========================================== c subroutine wrtfac(e,z1,nr,n,nvec,ni,io) c************************************************* c n - Dimension of sample covariance, nvec - Number of common factors c Eigenvalues of (er**-1/2)*cov*(er**-1/2) => e(n) c================================================= real*8 e(n),ch,z1 intent(in) e,z1,nr,n,nvec,ni,io c************************************************* c Output eigenvectors write(io,'(/7x,''Maximum Likelihood Factor Analysis Eigenvalues'') *') write(io,'(/11x,''Eigenvalue'',5x,''%'',9x,''% Total''/)') t=0. sum=0. do k=1,n sum=sum+sngl(e(k)) enddo do k=1,n s=100.*sngl(e(k))/sum t=t+s write(io,'(1x,i2,2x,g14.2,4x,2(f6.2,7x))') k,e(k),s,t enddo z2=dfloat(2*n+5)/6.d0+dfloat(2*nvec)/3.d0 nd=((n-nvec)*(n-nvec)-(n+nvec))/2 ch=z1*(dfloat(nr)-z2) write(io,'(/2x,''Chi Square Statistic = '',f12.1,2x,''Degrees of F *reedom = '',i4)') sngl(ch),nd write(io,'(7x,''Number of Iterations = '',i4)') ni c------- return end c c=========================================== c subroutine phi(cov,a,v,e,er,w,n,nvec,io) c************************************************* c nvec - Number of common factors, cov - Sample covariance matrix c er(n) - Specific variances c e(n) - Eigenvalues (er**-1/2)*cov*(er**-1/2), v - eigenvectors c================================================= implicit real*8(a-h,o-z) real*8 a(n,n),cov(n,n),v(n,n),e(n),er(n),w(n) intent(in) n,cov,er,io intent(out) a,e c************************************************* c a=(er**-1/2)*cov*(er**-1/2) do i=1,n do j=i,n a(j,i)=cov(j,i)/dsqrt(er(i)*er(j)) a(i,j)=a(j,i) enddo enddo c Eigenvalues & eigenvectors of a call eigen(a,n,n,e,w,v,jev) c a=(er**-1)-(er**-1/2)*(vv')*(er**-1/2) do i=1,n c Trap negative e's if(e(i)<=0.) then e(i)=1.d0 endif do j=i,n z=0.d0 do k=1,nvec z=z+v(i,k)*v(j,k) enddo a(j,i)=-z/dsqrt(er(i)*er(j)) a(i,j)=a(j,i) enddo a(i,i)=a(i,i)+1.d0/er(i) enddo c a**2 do i=1,n do j=i,n a(j,i)=a(j,i)*a(j,i) a(i,j)=a(j,i) enddo enddo c------- return end c c=========================================== c subroutine fd(cov,v,e,er,n,nvec,f) c******************************************** c 1st derivatives f(i)=del(S)/del(er(i)) c============================================ implicit real*8(a-h,o-z) real*8 cov(n,n),v(n,n),e(n),er(n),f(n) intent(in) cov,v,e,er,n,nvec intent(out) f c******************************************** do i=1,n z=0.d0 do j=1,nvec z=z+(e(j)-1.d0)*(v(i,j)*v(i,j)) enddo f(i)=(z+1.d0-cov(i,i)/er(i))/er(i) enddo c------- return end c c=========================================== c subroutine sd(v,e,er,n,nvec,f,a,d) c******************************************** c 2nd derivative matrix of S => d c============================================ implicit real*8(a-h,o-z) real*8 v(n,n),e(n),er(n),f(n),a(n,n),d(n,n) integer r intent(in) a,e,er,f,v,n,nvec intent(out) d c******************************************** nc=nvec+1 do i=1,n do j=i,n z=0.d0 do k=1,nvec y=0.d0 do r=nc,n y=y+2.d0*e(k)*(e(r)-1.d0)*v(i,r)*v(j,r)/(e(r)-e(k)) enddo z=z+y*v(i,k)*v(j,k) enddo d(j,i)=a(j,i)+z d(i,j)=d(j,i) enddo d(i,i)=d(i,i)-2.d0*f(i)/er(i) enddo c------- return end c c=========================================== c subroutine dludcmp(a,n,np,indx,d,ier,io) c************************************************* c Computes the LU decomposition of the matrix a c================================================= implicit real*8(a-h,o-z) parameter(nmax=500,tiny=1.0d-20) real*8 a(np,np),vv(nmax) real d integer indx(n),nmax intent(in) n,np,io intent(out) d,ier c************************************************* d=1. ier=0 do i=1,n aamax=0.d0 do j=1,n if(dabs(a(i,j))>aamax) then aamax=dabs(a(i,j)) endif enddo if(aamax==0.d0) then write(io,*) 'Singular matrix in ludcmp' ier=1 return endif vv(i)=1.d0/aamax enddo do j=1,n do i=1,j-1 sum=a(i,j) do k=1,i-1 sum=sum-a(i,k)*a(k,j) enddo a(i,j)=sum enddo aamax=0.d0 do i=j,n sum=a(i,j) do k=1,j-1 sum=sum-a(i,k)*a(k,j) enddo a(i,j)=sum dum=vv(i)*dabs(sum) if(dum>=aamax) then imax=i aamax=dum endif enddo if(j/=imax) then do k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum enddo d=-d vv(imax)=vv(j) endif indx(j)=imax if(dabs(a(j,j))0) then ir=ir+1 do i=1,mn z(i,ir)=va(i) enddo endif c------- c OLS fit - no intercept call regres(e,z,b,er,mn,ir,0,r2,0,ier,io) c Solve for parameters bs=1./b(1) x1b=0. if(cns==0) then x1b=-0.5/b(2) endif if(bs<=1.e-7) then write(io,*) 'Warning' write(io,*) 'bs = ',bs,' < 0. Absolute value is used' endif bs=sqrt(abs(bs)) x1b=x1b*bs c Estimate of candidate positions do k=1,m cr(k,1)=v(k)/bs enddo if(vi>0) then wv=b(ir) else wv=0. endif c------- return end c c=========================================== c subroutine param(b,cns,b1,b2,cosr,sinr,r,ro,xm,io) c******************************************** c Parameters from least squares fit of mean score c ro - rotation angle in radians, cosr - cosine(ro) c b1 & b2 - sigmas of x1 & x2, xm - mean ideal point c r - fixed rotation angle if cns=2, xm=0 if cns=1 c============================================ real*8 b(5),del,s1,s2 real xm(2),r integer cns intent(in) b,cns,r,io intent(out) b1,b2,cosr,sinr,ro,xm c******************************************** pi=3.141592 s1=b(2)-b(1) s2=b(1)+b(2) if(dabs(s1)>1.d-14) then del=datan2(b(3),s1)/2.d0 ro=sngl(del) del=dsin(2.d0*del) b1=sngl(s2-b(3)/del)/2.0 b2=sngl(s2+b(3)/del)/2.0 if(cns/=2) then c Rotation cosine & sin cosr=cos(sngl(del)) sinr=sin(sngl(del)) endif else if(cns/=2) then cosr=1.0/sqrt(2.0) sinr=cosr ro=pi/2.0 b1=sngl(s2-b(3))/2.0 b2=sngl(s2+b(3))/2.0 endif endif if(cns==2) then c Rotation set to angle r angle=pi*r/180.0 ro=angle cosr=cos(angle) sinr=sin(angle) endif c------- c Trap negative b's if(b1<=1.e-7) then write(io,*)'Warning' write(io,*)'b1 = ',b1,' < 0. b1 set to 1.e-4' b1=1.d-4 endif if(b2<=1.e-7) then write(io,*)'Warning' write(io,*)'b2 = ',b2,' < 0. b2 set to 1.e-4' b2=1.d-4 endif c Mean ideal points x1b=(b(4)*cosr-b(5)*sinr)/2. x2b=(-b(4)*sinr-b(5)*cosr)/2. b1=1./sqrt(b1) b2=1./sqrt(b2) if(cns==0) then xm(1)=x1b*b1 xm(2)=x2b*b2 endif if(cns==1) then c Mean=0 xm(1)=0. xm(2)=0. endif c------- return end c c=========================================== c subroutine v2(m,c,cr,z,d,e,b,er,vi,va,wv,cns,r,r2,b0,b1,b2,xc,xm,x *0,io) c******************************************** c Estimate candidate positions in two dimensions. (Note: array c is fc in call) c m - Candidates used (mu in call), c - 2 columns of rotated factors <= rotate c vi = 0 - No candidate characteristic values used c = 1 - Candidate characteristic weight estimated c va - Candidate characteristic values derived from survey data c cns = 0 All parameters estimated in OLS fit c = 1 Set the origin of the map => x0 from input xc c = 2 - Set rotation => cr c b0 - Rotation angle. If cns = 2 then b0=r. c e - adjusted dbar c Output: cr - Rotated candidate positions, wv - valence weight, r2 - r**2 c xm - mean ideal point, x0 - scaled input xc c============================================ integer m,vi,cns,sw real*8 e(m),b(m),d(m-1,m-1),z(m-1,6),er(m),del3,s1,s2 real c(m-1,2),cr(m,2),va(m),xc(2),xm(2),x0(2),b0,b1,b2,pi,r,an,r2, *wv intent(in) m,c,e,vi,va,cns,r,xc,io intent(out) b0,b1,b2,cr,r2,wv,xm,x0 c******************************************** b0=0.0 b1=0.0 b2=0.0 mn=m-1 ir=5 c Set variables for OLS fit do k=1,mn z(k,1)=dble(c(k,1))*dble(c(k,1)) z(k,2)=dble(c(k,2))*dble(c(k,2)) z(k,3)=dble(c(k,1))*dble(c(k,2)) z(k,4)=dble(c(k,1)) z(k,5)=dble(c(k,2)) enddo c Check for valience weight if(vi==1) then ir=6 endif do k=1,mn z(k,6)=va(k) enddo c------- pi=3.141592 c OLS fit int=0 sw=1 call regres(e,z,b,er,mn,ir,int,r2,sw,ier,io) c Solve for parameters s1=b(2)-b(1) s2=b(1)+b(2) if(dabs(s1)>1.d-7) then del3=datan2(b(3),s1)/2.d0 b0=180.0*sngl(del3)/pi if(cns<=1) then cos3=cos(sngl(del3)) sin3=sin(sngl(del3)) endif del3=dsin(2.d0*del3) b1=sngl(b(1)+b(2)-b(3)/del3)/2.0 b2=sngl(b(1)+b(2)+b(3)/del3)/2.0 c------- c Trap negative b's if(b1<=1.e-7) then write(io,'(/2x,''Warning b1 = '',f7.2,'' < 0. b1 set to 1.'')')b1 b1=1.0 endif if(b2<=1.e-7) then write(io,'(/2x,''Warning b2 = '',f7.2,'' < 0. b2 set to 1.'')')b2 b2=1.0 endif b1=1.0/sqrt(b1) b2=1.0/sqrt(b2) c------- x0(1)=0.0 x0(2)=0.0 if(cns==1) then c Constrain origin => xc from input xc x1b=xc(1) x2b=xc(2) x0(1)=x1b*b1 x0(2)=-x2b*b2 endif c Mean ideal point if(cns/=1) then x1b=(b(4)*cos3-b(5)*sin3)/2.0 x2b=(b(4)*sin3+b(5)*cos3)/2.0 xm(1)=x1b*b1 xm(2)=-x2b*b2 endif c------- if(cns==2) then c Contrain rotation b0=r an=r*(pi/180.0) cos3=cos(an) sin3=sin(an) endif else c If b(1)=b(2) if(s1==0.d0) then cos3=1.0/sqrt(2.0) sin3=cos3 endif b0=45.0 b1=sngl(s2-b(3))/2.0 b2=sngl(s2+b(3))/2.0 endif c------- c Candidate positions do k=1,mn s1=c(k,1)*cos3-c(k,2)*sin3 s2=c(k,1)*sin3+c(k,2)*cos3 cr(k,1)=s1/b1-x0(1) cr(k,2)=s2/b2-x0(2) enddo c------- if(vi==1) then wv=b(ir) else wv=0. endif c Candidate m cr(m,1)=-x0(1) cr(m,2)=-x0(2) c------- return end c c=========================================== c subroutine v1(m,c,cr,z,d,e,b,er,vi,va,wv,cns,b1,r2,io) c******************************************** c Estimate candidate positions in one dimension.(c is yc in call) c m - Number of candidates used (mu in call), c - Rotated factors from rotate c vi= 0 - No candidate characteristic values used c = 1 - Candidate characteristic weight estimated c va - Candidate characteristic values derived from survey data c cns(1)= All parameters estimated in OLS fit c = 1 Set mean ideal point c Output: cr - Rotated candidate positions, wv - valence weight, r2 - r**2 c================================================= integer vi,cns(2) real*8 e(m),b(m),d(m-1,m-1),z(m-1,6),er(m) real c(m),cr(m,1),va(m),b1,r2,wv intent(in) m,c,e,vi,va,cns,io intent(out) cr,b1,r2,wv c************************************************* mn=m-1 ir=2 c Check for constraints on mean if(cns(1)==1) then ir=1 endif c Set variables for OLS fit do k=1,mn z(k,1)=dble(c(k))*dble(c(k)) z(k,2)=dble(c(k)) z(k,3)=1.0 enddo c Check for valience weight if(vi==1) then ir=ir+1 do i=1,mn z(i,ir)=va(i) enddo endif c------- c OLS fit call regres(e,z,b,er,mn,ir,1,r2,2,ier,io) c Solve for parameters b1=1.0/b(1) x1b=0.0 if(cns(1)==0) then x1b=-0.5/b(2) endif if(b1<=1.e-7) then write(io,*)'Warning' write(io,*)'b1 = ',b1,' < 0. Absolute value is used' endif b1=sqrt(abs(b1)) x1b=x1b*b1 c Estimate of candidate positions do k=1,mn cr(k,1)=c(k)/b1-x1b enddo if(vi==1) then wv=b(ir) else wv=0.0 endif cr(m,1)=-x1b c------- return end c c=========================================== c subroutine regres(y,x,b,c,n,p,int,r2,sw,ier,io) c************************************************* c Input: x - matrix of independent variables, y - dependent variables c n - No. of observations, p - No. of independent variables c int - Fit a mean if int = 1, sw=2 - output stats, sw=1 & 2 - header c Output: b - Regression coefficients, c - standard errors, r2 c================================================= implicit real*8(a-h,o-z) allocatable::a,s,w integer p,n,int,sw real r2 real*8 y(n),x(n,p),b(p),c(p),a(:,:),w(:),s(:) intent(in) y,x,n,p,int,sw,io intent(out) b,c,r2,ier c************************************************* pp=p-1 allocate(a(0:pp,0:pp),s(0:p*p-1),w(0:pp),stat=ibad) if(ibad/=0) then write(io,*)'Unable to allocate work space for s in regres' return endif sy=0.d0 do k=0,pp w(k)=0.d0 enddo if(int==1) then c Intercept rn=dfloat(n) do i=1,n sy=sy+y(i) do k=1,p w(k-1)=w(k-1)+x(i,k) enddo enddo c Mean y sy=sy/rn do k=0,pp c Mean of x(i,k) w(k)=w(k)/rn enddo endif c------- do j=1,p jj=j-1 do k=j,p kk=k-1 a(jj,kk)=0.d0 enddo enddo c Covariance matrix x*x do i=1,n do j=1,p jj=j-1 dij=x(i,j)-w(jj) a(jj,jj)=a(jj,jj)+dij*dij do k=j+1,p kk=k-1 dij=x(i,j)-w(jj) dik=x(i,k)-w(jj) a(jj,kk)=a(jj,kk)+dij*dik a(kk,jj)=a(jj,kk) enddo enddo enddo c Cov(x,y) do k=1,p st=0.d0 do i=1,n st=st+(x(i,k)-w(k-1))*(y(i)-sy) enddo c(k)=st enddo c Sum of squared y-mean(y) ssy=0.d0 do i=1,n ssy=ssy+(y(i)-sy)*(y(i)-sy) enddo c------- c a=>s do i=1,p ii=i-1 do j=1,p jj=j-1 st=a(jj,ii) s((i-1)*p+j-1)=st enddo enddo call dls(s,b,c,w,p,1,ier,io) if(ier==1) then write(io,*)'Covariance matrix in regres is not positive definite' return endif ca=0.d0 if(int==1) then c Mean y do k=1,p ca=ca+w(k-1)*b(k) enddo endif c Intercept cr cr=sy-ca r=0.d0 c Estimate y's do i=1,n st=0.d0 do j=1,p st=st+x(i,j)*b(j) enddo c Residuals res=y(i)-sy-st r=r+res*res enddo c------- sdy=dsqrt(ssy/dfloat(n-1)) r2=sngl(1.d0-r/ssy) r2=r2-float((p)/(n-p-1))*(1-r2) se=dsqrt(r/dfloat(n-p)) if(sw==1.or.sw==2) then c Output statistics write(io,'(2x,65(''=''))') write(io,'(12x,''Statistics of OLS Fit''/)') write(io,'(2x,''Standard Deviation of Y ='',f15.3)')sdy write(io,'(2x,''R Square = '',f7.4,2x,''Standard Error = '',f15.4 *)')r2,se write(io,'(2x,65(''+''))') endif if(int==1) then write(io,'(7x,''Intercept = '',g10.3/)')cr endif c Standard errors of OLS estimates do k=1,p c(k)=se*dsqrt(s((k-1)*p+k-1)) enddo if(sw==2) then write(io,'(2x,65(''=''))') write(io,'(12x,''Statistics of OLS Fit''/)') write(io,'(/7x,''OLS Estimates'',6x,''Standard Error''/)') write(io,'(3x,i2,2x,g10.3,10x,g10.3)')(k,b(k),c(k),k=1,p) write(io,'(2x,65(''+''))') endif c------- deallocate(a,s,w) return end c c=========================================== c subroutine eigen(cov,np,n,d,e,z,jev) c******************************************** c Finds eigenvalues and eigenvectors of a symmetric real matrix COV c whose order is n. np is the row dimension of cov set in the c program dimension (or parameter) statement c Output: c d contains the eigenvalues in descending order. If jev = j, then c eigenvalues are correct but unordered for indices 1,2,..,jev-1. c jth column of z is the j-th orthonormal eigenvector. c jev = zero for normal return, c = j if the j-th eigenvalue has not been determined after c 30 iterations. c============================================ real*8 cov(np,n),d(n),e(n),z(np,n) intent(in) cov,np,n intent(out) d,z,jev c******************************************** if(n<1) then write(6,*) 'Matrix order n < 2 in subroutine eigen' return elseif(n>np) then write(6,*) 'Matrix order n > limit',np,'in subroutine eigen' return endif c------- c Call tridiagonalizing subroutine call tridag(cov,d,e,z,n,np) c Call main eigenvalue program call evs(d,e,z,n,np,jev) c-------- return end c c============================================ c subroutine tridag(a,d,e,z,n,np) c******************************************** c This subroutine reduces a real symmetric matrix to a symmetric c tridiagonal matrix using orthogonal similarity transformations. c Input: c n is the order of the matrix; c a contains the real*8 symmetric input matrix; c Output: c d contains the diagonal elements of the tridiagonal matrix; c e contains the subdiagonal elements of the tridiagonal c matrix in its last n-1 positions. e(1) is set to zero; c z contains the orthogonal transformation matrix c produced in the reduction; c============================================ integer i,j,k,l,n,ii,jp1 real*8 a(np,n),d(n),e(n),z(np,n),f,g,h,hh,scale c******************************************** do i=1, n d(i)=0.d0 do j = 1, i z(i,j)=a(i,j) enddo enddo do 30 ii = 2, n i = n + 2 - ii l = i - 1 h = 0.d0 scale = 0.d0 if (l < 2) go to 17 c :::::::::: scale row :::::::::: do 12 k = 1, l 12 scale = scale + dabs(z(i,k)) c if (scale /= 0.d0) go to 14 17 e(i) = z(i,l) go to 290 c 14 do 15 k = 1, l z(i,k) = z(i,k) / scale h = h + z(i,k) * z(i,k) 15 continue c f = z(i,l) g = -dsign(dsqrt(h),f) e(i) = scale * g h = h - f * g z(i,l) = f - g f = 0.d0 c do 24 j = 1, l z(j,i) = z(i,j) / h g = 0.d0 c :::::::::: form element of a*u :::::::::: do 180 k = 1, j 180 g = g + z(j,k) * z(i,k) c jp1 = j + 1 if (l < jp1) go to 220 c do 20 k = jp1, l 20 g = g + z(k,j) * z(i,k) c :::::::::: form element of p :::::::::: 220 e(j) = g / h f = f + e(j) * z(i,j) 24 continue c hh = f / (h + h) c :::::::::: form reduced a :::::::::: do 260 j = 1, l f = z(i,j) g = e(j) - hh * f e(j) = g c do 260 k = 1, j z(j,k) = z(j,k) - f * e(k) - g * z(i,k) 260 continue c 290 d(i) = h 30 continue c :::::::::: accumulation of transformation matrices :::::::::: do 500 i = 1, n l = i - 1 if (d(i) == 0.d0) go to 380 c do 360 j = 1, l g = 0.d0 c do 340 k = 1, l 340 g = g + z(i,k) * z(k,j) c do 360 k = 1, l z(k,j) = z(k,j) - g * z(k,i) 360 continue c 380 d(i) = z(i,i) z(i,i) = 1.d0 if (l < 1) go to 500 c do 400 j = 1, l z(i,j) = 0.d0 z(j,i) = 0.d0 400 continue c 500 continue c------- return end c c============================================ c subroutine evs(d,e,z,n,np,jev) c******************************************** c Input: c n is the order of the matrix; c d contains the diagonal elements of the input matrix; c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary; c z contains the transformation matrix produced in the c reduction by TRIDAG if called. If the eigenvectors c of the tridiagonal matrix are desired, z must contain c the identity matrix. c Output: c d contains the eigenvalues in descending order. If an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,jev-1; c e has been destroyed; c z contains orthonormal eigenvectors of the symmetric c tridiagonal (or full) matrix. If an error exit is made, c z contains the eigenvectors associated with the stored ev's; c jev = zero for a normal return c = j if the j-th eigenvalue has not been determined c after 30 iterations. c============================================ integer i,j,k,l,m,n,ii,mml,jev real*8 d(n),e(n),z(np,n),b,c,f,g,p,r,s,machep c--------------------------------- c machep is a machine dependent parameter specifying c the relative precision of floating point arithmetic. c machep = 16.0d0**(-13) for long form arithmetic on s360 data machep/16.d-40/ c******************************************** jev = 0 do 10 i = 2, n 10 e(i-1) = e(i) c e(n) = 0.d0 c do 24 l = 1, n j = 0 c :::::::::: look for small sub-diagonal element :::::::::: 1 do 11 m = l, n if (m == n) go to 12 if (dabs(e(m))<=machep * (dabs(d(m))+dabs(d(m+1)))) goto 12 c 11 continue c 12 p = d(l) if (m == l) goto 24 c set error -- no convergence to an eigenvalue after 30 iterations if (j == 30) then jev=1 return endif j = j + 1 c :::::::::: form shift :::::::::: g = (d(l+1) - p) / (2.d0 * e(l)) r = dsqrt(g*g+1.d0) g = d(m) - p + e(l) / (g + dsign(r,g)) s = 1.d0 c = 1.d0 p = 0.d0 mml = m - l c :::::::::: for i=m-1 step -1 until l do -- :::::::::: do 20 ii = 1, mml i = m - ii f = s * e(i) b = c * e(i) if (dabs(f) < dabs(g)) goto 15 c = g / f r = dsqrt(c*c+1.d0) e(i+1) = f * r s = 1.d0 / r c = c * s goto 16 15 s = f / g r = dsqrt(s*s+1.d0) e(i+1) = g * r c = 1.d0 / r s = s * c 16 g = d(i+1) - p r = (d(i) - g) * s + 2.d0 * c * b p = s * r d(i+1) = g + p g = c * r - b c :::::::::: form vector :::::::::: do 18 k = 1, n f = z(k,i+1) z(k,i+1) = s * z(k,i) + c * f z(k,i) = c * z(k,i) - s * f 18 continue c 20 continue c d(l) = d(l) - p e(l) = g e(m) = 0.d0 go to 1 24 continue c :::::::::: order eigenvalues and eigenvectors :::::::::: do 30 ii = 2, n i = ii - 1 k = i p = d(i) c do 26 j = ii, n if (d(j) > p) then k = j p = d(j) endif 26 continue c if (k /= i) then d(k) = d(i) d(i) = p c do 28 j = 1, n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 28 continue endif c 30 continue c------- return end c c============================================ c subroutine dls(a,x,b,d,n,inv,ier,io) c******************************************** c Solves the n linear equations Ax=b c Input: a - positive definite matrix c Output: x - solution, a - inverse of a matrix if inv = 1 c============================================ real*8 a(n,n),b(n),d(n),x(n),sum intent(in) b,n,inv,io intent(out) x,ier intent(in out) a c******************************************** ier=0 do i=1,n do j=i,n sum=a(i,j) do k=i-1,1,-1 sum=sum-a(i,k)*a(j,k) enddo if(i==j) then if(sum<=0.d0) then ier=1 c write(io,*) 'a is not positive definite in DLS' return endif d(i)=dsqrt(sum) else a(j,i)=sum/d(i) endif enddo enddo c------- do i=1,n sum=b(i) do k=i-1,1,-1 sum=sum-a(i,k)*x(k) enddo x(i)=sum/d(i) enddo do i=n,1,-1 sum=x(i) do k=i+1,n sum=sum-a(k,i)*x(k) enddo x(i)=sum/d(i) enddo if(inv/=1) then return endif c------- c Lower diagonal inverse of L in a do i=1,n a(i,i)=1.d0/d(i) do j=i+1,n sum=0.d0 do k=i,j-1 sum=sum-a(j,k)*a(k,i) enddo a(j,i)=sum/d(j) enddo enddo c------- c Inverse of LL' do i=1,n do j=i+1,n sum=0.d0 do k=j,n sum=sum+a(k,i)*a(k,j) enddo a(i,j)=sum enddo c-------- c Diagonal sum=0.d0 do k=i,n sum=sum+a(k,i)*a(k,i) enddo a(i,i)=sum enddo do i=1,n do j=i+1,n a(j,i)=a(i,j) enddo enddo c------- return end c c=========================================== c subroutine ideal2(n,nu,m,x,uv,c,dbar,va,wv,a,e,w,y,vd,dav,pow,ref, *sc,io) c******************************************** c Fit ideal points in a two dimensional space. Version 5-7-2008 c x - Input: Transformed rearranged score columns (c in call) c Best score is 0 c vd - respondent id's, dbar - mean candidate score, m - no. candidates used c uv - flag for candidate whose votes are to be estimated c Output: Ideal points. x(i,3) = kt if one or more candidates have a 0 score, c x(i,3) = 11 if the line segment fit is used, x(i,3) = 21 if the least c squares fit is used c ni - number used, a - Sum of voters nearest to each candidate c c - Matrix of candidate positions (cr in call) c xmn =>y - mean of estimated ideal points c============================================ real x(nu,m),c(m,2),va(m),y(7),wv,sc real*8 a(m),dbar(m),e(m) integer uv(m),w(m),pow character*5 vd(nu) character*2 ref character dav intent(in) n,nu,m,c,dbar,va,wv,vd,dav,pow,ref,sc,uv,io intent(out) a,y intent(in out) x c************************************************* mn=m-1 do k=1,5 y(k)=0.0 enddo c X'X of candidate points do k=1,mn smx=c(k,1)-c(m,1) smn=c(k,2)-c(m,2) y(1)=y(1)+smx*smx y(2)=y(2)+smn*smn y(3)=y(3)+smx*smn enddo c Determinant*2 rm=(y(1)*y(2)-y(3)*y(3))*2. c------- c Initialize seed call random_seed xm1=0.0 xm2=0.0 do i=1,nu c Add back last column do j=1,mn e(j)=dble(x(i,j)+x(i,m)) enddo e(m)=dble(x(i,m)) c Find minimum & maximum of the ith row transformed scores (smn,smx) smn=sngl(e(1)) smx=smn kt=0 do j=1,m st=sngl(e(j)) if(stsmx) then smx=st endif c------- c Find j's for zero transformed scores if(st==0.0) then c Count zero scores kt=kt+1 w(kt)=j endif enddo c Set ideal point at the candidate with a 0. If 11) then x1=0. x2=0. do k=1,kt x1=x1+c(w(k),1) x2=x2+c(w(k),2) enddo x(i,1)=x1/float(kt)+sc*(y(6)-0.5) x(i,2)=x2/float(kt)+sc*(y(7)-0.5) x(i,3)=float(min(kt,10)) endif c-------- c Convex combination on line kmn=0 kmx=0 if(kt==0) then do j=1,m as=sngl(e(j)) if(abs(as-smn)<1.e-7) then c j has the lowest score at the min value of transformed grades kmn=kmn+1 w(1)=j endif if(abs(as-smx)<1.e-7.and.as>0.0) then c j has the highest score but not max value kmx=kmx+1 w(2)=j endif enddo c------- if(kmn==1.and.kmx==1) then if(pow==2) then a1=sqrt(sngl(e(w(1)))) a2=sqrt(sngl(e(w(2)))) else a1=sngl(e(w(1))) a2=sngl(e(w(2))) endif c Ratio of inverse grades st=-(a2/a1) dc=(c(w(1),1)-c(w(2),1))**2+(c(w(1),2)-c(w(2),2))**2 ar=st/sqrt(dc) c Prorate on cord between thetas call random_number(y(6:7)) x(i,1)=ar*c(w(2),1)+(1.-ar)*c(w(1),1)+sc*(y(6)-.5) x(i,2)=ar*c(w(2),2)+(1.-ar)*c(w(1),2)+sc*(y(7)-.5) x(i,3)=11. else c Least squares fit using the adjusted scores y(4)=0.0 y(5)=0.0 do k=1,mn e(k)=dble(x(i,k)-dbar(k)) y(4)=y(4)+(c(k,1)-c(m,1))*sngl(e(k)) y(5)=y(5)+(c(k,2)-c(m,2))*sngl(e(k)) enddo c Ideal points x(i,1)=(y(3)*y(5)-y(2)*y(4))/rm+c(m,1) x(i,2)=(y(3)*y(4)-y(1)*y(5))/rm+c(m,2) x(i,3)=21. endif endif enddo c------- write(io,'(2x,70(''*''))') do k=1,m a(k)=0.d0 enddo c Find minimum distance from the kth voter to the candidates do k=1,nu smn=1.e7 knc=0 do nc=1,m if(uv(nc)==1) then c Distance from x(k) to c(nc) distc=(c(nc,1)-x(k,1))**2+(c(nc,2)-x(k,2))**2+wv*va(nc) if(distc<0.) then c Dist=0 dist=0. endif if(distc>=0.) then dist=sqrt(distc) endif if(distsmx) then smx=st endif c------- c Find j's for zero scores if(st==0.0) then kt=kt+1 w(kt)=j endif enddo c End of j loop call random_number(y) c Set ideal point at the candidate with a 0. If 11) then x1=0.0 do k=1,kt x1=x1+c(w(k)) enddo x(i,1)=x1/float(kt)+sc*y x(i,2)=float(min(kt,10)) endif if(kt==0) then do j=1,m st=x(i,j) if(abs(st-smn)<1.e-7) then w(1)=j endif if(abs(st-smx)<1.e-7) then w(2)=j endif enddo c------- c Ideal point from quadratic model d=c(w(1))-c(w(2)) st=(x(i,w(1))-x(i,w(2)))/(2.0*d*d) x(i,1)=(0.5-st)*c(w(1))+(0.5+st)*c(w(2))+sc*y x(i,2)=11. endif enddo c------- return end c c=========================================== c subroutine readmap(id,n,m,kis,mpflg,fmt,ii,io) c******************************************** c Read id, n - respondents, m - candidates, kis - issue flags columns, c mpflg - map flags=0/1, fmt - format c============================================ character*120 buf character*80 id character*50 charead,fmt character*20 par character*2 delim intent(in) ii,io intent(out) id,n,m,kis,mpflg,fmt c**************************** rewind(ii) c Read file header read(ii,'(a80)')id read(ii,'(a120)')buf c Sample size (respondents) delim(1:2)='= ' par='ndents' n=numread(par,buf,delim,ir,ier,io) if(ier>0) then write(io,*)'Error in reading the no. of respondents' write(io,*)buf write(6,*)'Error in reading the no. of respondents' write(6,*)buf write(9,*)'Error in reading the no. of respondents' write(9,*)buf stop endif c------- c Candidates par='dates' m=numread(par,buf,delim,ir,ier,io) if(ier>0) then write(io,*)'Error in reading the number of candidates' write(io,*)buf write(6,*)'Error in reading the number of candidates' write(6,*)buf write(9,*)'Error in reading the number of candidates' write(9,*)buf stop endif c------- c Issue means flag columns kis=0 par='olumns' kis=numread(par,buf,delim,ir,ier,io) if(ier>0) then write(io,*)'Error in reading the no. of issue flag columns' write(io,*)buf write(6,*)'Error in reading the no. of issue flag columns' write(6,*)buf write(9,*)'Error in reading the no. of issue flag columns' write(9,*)buf return endif c------- c Map respondent column mpflg=0 par='ion flags' mpflg=numread(par,buf,delim,ir,ier,io) if(ier>0.or.(mpflg<0.or.mpflg>1)) then write(io,*)'Error in reading respondent selection flags switch' write(io,*)'map column no = ',mpflg write(6,*)'Error in reading respondent selection flags switch' write(6,*) 'map column no = ',mpflg write(9,*)'Error in reading respondent selection flags switch' write(9,*) 'map column no = ',mpflg return endif c------- c Format par='rmat' fmt=charead(par,buf,delim,ir,ier,3) if(ier>0) then write(io,*)'Error in reading the format from the data file' write(io,*)buf write(6,*)'Error in reading the format from the data file' write(6,*)buf write(9,*)'Error in reading the format from the data file' write(9,*)buf stop endif c------- return end c c============================================= c subroutine flagmed(x,rflg,m,kis,n,nu,flag,ng,check,h,med,i2,io) c******************************************** c Means of x's using flag & write them to io & ii if i2>0 c ng - no. of flag groups for kth flag line c============================================ allocatable::kt,nar real x(nu,2),h(nu,kis*12),med(kis,12),fi(1),fo(1) integer flag(kis,6),rflg(n,kis),ng(kis),kt(:),nar(:),t character*30 check(m+kis) character*20 name character acm,acn intent(in) x,check,m,kis,n,nu,flag,ng,i2,io intent(out) med c******************************************** allocate(kt(0:kis*6-1),nar(0:nu-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate work space for nar in flagmed' stop endif do k=0,kis*6-1 kt(k)=0 enddo do t=1,nu do kf=1,kis if(ng(kf)==1) then c One group na=(kf-1)*6 if((rflg(t,kf)==flag(kf,1).or.rflg(t,kf)==flag(kf,2)))then c Keep if flags check kt(na)=kt(na)+1 do k=1,2 kc=(kf-1)*2+k h(kt(na),kc)=dble(x(t,k)) enddo endif endif if(ng(kf)==2) then na=(kf-1)*6+1 c Two groups if((rflg(t,kf)==flag(kf,1).or.rflg(t,kf)==flag(kf,2))) then c Keep if flags check for group 1 kt(na)=kt(na)+1 do k=1,2 kc=kis*2+(kf-1)*2+k h(kt(na),kc)=dble(x(t,k)) enddo endif if((rflg(t,kf)==flag(kf,3).or.rflg(t,kf)==flag(kf,4))) then c Keep if flags check for group 2 ns=na+1 kt(ns)=kt(ns)+1 do k=1,2 kc=kis*4+(kf-1)*2+k h(kt(ns),kc)=dble(x(t,k)) enddo endif endif if(ng(kf)==3) then c Three groups na=(kf-1)*6+3 if((rflg(t,kf)==flag(kf,1).or.rflg(t,kf)==flag(kf,2))) then c Keep if flags check for group 1 kt(na)=kt(na)+1 do k=1,2 kc=kis*6+(kf-1)*2+k h(kt(na),kc)=dble(x(t,k)) enddo endif if((rflg(t,kf)==flag(kf,3).or.rflg(t,kf)==flag(kf,4))) then c Keep if flags check for group 2 ns=na+1 kt(ns)=kt(ns)+1 do k=1,2 kc=kis*8+(kf-1)*2+k h(kt(ns),kc)=dble(x(t,k)) enddo endif if((rflg(t,kf)==flag(kf,5).or.rflg(t,kf)==flag(kf,6))) then c Keep if flags check for group 3 ns=na+2 kt(ns)=kt(ns)+1 do k=1,2 kc=kis*10+(kf-1)*2+k h(kt(ns),kc)=dble(x(t,k)) enddo endif endif enddo c End of flag check loop enddo c End of respondent loop c------- if(i2>0) then write(i2,'(17x,''Issue or Type Medians''/)') endif ier=0 fi(1)=0.5 do kf=1,kis if(ng(kf)==1) then na=(kf-1)*6 nm=kt(na) if(nm>9) then do k=1,2 kc=(kf-1)*2+k call sort(nar,h(1,kc),nm,1,fi,fo,'B',ier) if(ier==1) then write(6,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf write(9,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf endif med(kf,k)=fo(1) enddo else med(kf,1)=-1. med(kf,2)=-1. endif endif if(ng(kf)==2) then na=(kf-1)*6+1 nm=kt(na) if(nm>9) then do k=1,2 kc=kis*2+(kf-1)*2+k call sort(nar,h(1,kc),nm,1,fi,fo,'B',ier) if(ier==1) then write(6,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf write(9,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf endif med(kf,2+k)=fo(1) enddo else med(kf,3)=-1. med(kf,4)=-1. endif nm=kt(na+1) if(nm>9) then do k=1,2 kc=kis*4+(kf-1)*2+k call sort(nar,h(1,kc),nm,1,fi,fo,'B',ier) if(ier==1) then write(6,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf write(9,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf endif med(kf,4+k)=fo(1) enddo else med(kf,5)=-1. med(kf,6)=-1. endif endif if(ng(kf)==3) then na=(kf-1)*6+3 nm=kt(na) if(nm>9) then do k=1,2 kc=kis*6+(kf-1)*2+k call sort(nar,h(1,kc),nm,1,fi,fo,'B',ier) if(ier==1) then write(6,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf write(9,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf endif med(kf,6+k)=fo(1) enddo else med(kf,7)=-1. med(kf,8)=-1. endif nm=kt(na+1) if(nm>9) then do k=1,2 kc=kis*8+(kf-1)*2+k call sort(nar,h(1,kc),nm,1,fi,fo,'B',ier) if(ier==1) then write(6,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf write(9,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf endif med(kf,8+k)=fo(1) enddo else med(kf,9)=-1. med(kf,10)=-1. endif nm=kt(na+2) if(nm>9) then do k=1,2 kc=kis*10+(kf-1)*2+k call sort(nar,h(1,kc),nm,1,fi,fo,'B',ier) if(ier==1) then write(6,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf write(9,'(2x,''Length '',i3,'' sort array < 1 - line '',i2)') *nm,kf endif med(kf,10+k)=fo(1) enddo else med(kf,11)=-1. med(kf,12)=-1. endif endif c------- c Output issue ideal point means j=m+kf ia=index(check(j),' ') if(ia>18) then ia=18 endif name(:ia-1)=check(j)(:ia-1) do ks=ia+1,20 name(ks:ks)=' ' enddo c------- if(ng(kf)==1) then acm=achar(flag(kf,1)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,2)+48) name(ia+1:ia+1)=acn endif if(med(kf,1)/=-1.) then write(io,'(a20,1x,:,4(f12.3,1x))')name,(med(kf,k),k=1,2) endif endif if(ng(kf)==2) then acm=achar(flag(kf,1)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,2)+48) name(ia+1:ia+1)=acn endif if(med(kf,3)/=-1.) then write(io,'(a20,1x,:,4(f12.3,1x))')name,(med(kf,2+k),k=1,2) endif acm=achar(flag(kf,3)+48) name(ia:ia)=acm if(flag(kf,4)/=-1) then acn=achar(flag(kf,4)+48) name(ia+1:ia+1)=acn endif if(med(kf,5)/=-1.) then write(io,'(a20,1x,:,4(f12.3,1x))')name,(med(kf,4+k),k=1,2) endif endif if(ng(kf)==3) then acm=achar(flag(kf,1)+48) name(ia:ia)=acm if(flag(kf,1)/=-1) then acn=achar(flag(kf,2)+48) name(ia+1:ia+1)=acn endif if(med(kf,7)/=-1.) then write(io,'(a20,1x,:,4(f12.3,1x))')name,(med(kf,6+k),k=1,2) endif acm=achar(flag(kf,3)+48) name(ia:ia)=acm if(flag(kf,4)/=-1) then acn=achar(flag(kf,4)+48) name(ia+1:ia+1)=acn endif if(med(kf,9)/=-1.) then write(io,'(a20,1x,:,4(f12.3,1x))')name,(med(kf,8+k),k=1,2) endif acm=achar(flag(kf,5)+48) name(ia:ia)=acm if(flag(kf,4)/=-1) then acn=achar(flag(kf,6)+48) name(ia+1:ia+1)=acn endif if(med(kf,11)/=-1.) then write(io,'(a20,1x,:,4(f12.3,1x))')name,(med(kf,10+k),k=1,2) endif endif c------- if(i2>0) then mp=2 if(ng(kf)==1) then if(flag(kf,2)==-1) then mp=1 endif acm=achar(flag(kf,1)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,2)+48) name(ia+1:ia+1)=acn endif if(med(kf,1)/=-1.) then nm=kt(na) write(i2,'(a20,1x,2(f10.2,1x),2x,i5,'' voters for flag(s)= '', *:,2(i1,1x))')name,(med(kf,k),k=1,2),nm,(flag(kf,t),t=1,mp) endif endif c------- if(ng(kf)==2) then if(flag(kf,2)==-1) then mp=1 endif acm=achar(flag(kf,1)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,2)+48) name(ia+1:ia+1)=acn endif if(med(kf,3)/=-1.) then nm=kt((kf-1)*6+1) write(i2,'(a20,1x,2(f10.2,1x),2x,i5,'' voters for flag(s)= '', *:,2(i1,1x))')name,(med(kf,2+k),k=1,2),nm,(flag(kf,t),t=1,mp) endif if(flag(k,4)==-1) then mp=1 endif acm=achar(flag(kf,3)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,4)+48) name(ia+1:ia+1)=acn endif if(med(kf,5)/=-1.) then nm=kt((kf-1)*6+2) write(i2,'(a20,1x,2(f10.2,1x),2x,i5,'' voters for flag(s)= '', *:,2(i1,1x))')name,(med(kf,4+k),k=1,2),nm,(flag(kf,t+2),t=1,mp) endif endif c------- if(ng(k)==3) then if(flag(k,2)==-1) then mp=1 endif acm=achar(flag(kf,1)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,2)+48) name(ia+1:ia+1)=acn endif if(med(kf,7)/=-1.) then nm=kt((kf-1)*6+3) write(i2,'(a20,1x,2(f10.2,1x),2x,i5,'' voters for flag(s)= '', *:,2(i1,1x))')name,(med(kf,6+k),k=1,2),nm,(flag(kf,t),t=1,mp) endif if(flag(k,4)==-1) then mp=1 endif acm=achar(flag(kf,3)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,4)+48) name(ia+1:ia+1)=acn endif if(med(kf,9)/=-1.) then nm=kt((kf-1)*6+4) write(i2,'(a20,1x,2(f10.2,1x),2x,i5,'' voters for flag(s)= '', *:,2(i1,1x))')name,(med(kf,8+k),k=1,2),nm,(flag(kf,t+2),t=1,mp) endif if(flag(k,6)==-1) then mp=1 endif acm=achar(flag(kf,5)+48) name(ia:ia)=acm if(flag(kf,2)/=-1) then acn=achar(flag(kf,6)+48) name(ia+1:ia+1)=acn endif name(ia:ia)=acm if(med(kf,11)/=-1.) then nm=kt((kf-1)*6+5) write(i2,'(a20,1x,2(f10.2,1x),2x,i5,'' voters for flag(s)= '', *:,2(i1,1x))')name,(med(kf,10+k),k=1,2),nm,(flag(kf,t+4),t=1,mp) endif endif endif enddo c------- return end c c============================================= c subroutine lineread(buf,names,mu,io) c******************************************** c Read names in buf c============================================ character*700 buf character*30 names(mu) intent(in) buf,mu,io intent(out) names c******************************************** c Check for character in first space nc=ichar(buf(1:1)) if(nc/=32) then write(io,'(17x,''ERROR''/)') write(io,'(7x,''Leave at least one space before the first name in * the list of '',i2,'' candidate names'',/,6x,''above the first dat *a line'')')mu write(io,'(/(a80))')buf write(6,'(17x,''ERROR''/)') write(6,'(7x,''Leave at least one space before the first name in *the list of '',i2,'' candidate names'',/,6x,''above the first data * line'')')mu write(6,'(/(a80))')buf write(9,'(17x,''ERROR''/)') write(9,'(7x,''Leave at least one space before the first name in *the list of '',i2,'' candidate names'',/,6x,''above the first data * line'')')mu write(9,'(/(a80))')buf return endif npt=1 do nkt=1,mu c Search for first non blank kb=npt dowhile(ichar(buf(kb:kb))==32) kb=kb+1 enddo c Search for first blank after nkt name nc=kb dowhile(ichar(buf(nc:nc))/=32) nc=nc+1 enddo names(nkt)=buf(kb:nc-1) npt=nc enddo c------- return end c c=========================================== c subroutine new(can,cand,a,c,cov,cr,d,dbar,e,er,f,h,s,uc,uv,v,va,vs *,vd,w,x,flgd,m,mu,mn,n,nu,pow,dim,iw,psub,cns,dav,vi,fmo,ref,sc,r, *smis,nstr,del,smax,smin,vmax,xc,med,check,rflg,flag,kis,nfld,ng,nf *g,io) c******************************************** c uc - candidate selection 0/1, uv - vote count 0/1, va - valences c cns=1 - constrain mean, cns=2, r - rotation c============================================ real x(n,mu),c(nu,mu),cr(mu,2),va(m),vs(m),med(kis,12),y(7),bmn(4) *,br(2),bs1(2),bs2(2),bv(2),xm(2),xc(2),h(nu,kis*12),del real*8 a(mu*mu),cov(mn,mn),d(mn,mn),dbar(mu),e(mu),er(mu),f(mn),s( *mn*mn),v(mn,mn),w(mn*mn),z1 integer iw(mu),uc(m),uv(m),flag(kis,6),rflg(n,kis),ng(kis),flgd(10 *),vi,pow,dim,psub,cns,t character*50 fmo character*30 can(m),cand(m),check(m+kis) character*5 vd(n) character*2 ref character dav intent(in) can,cand,cov,flgd,uc,uv,va,m,mu,mn,n,nu,pow,dim,cns,vi, *fmo,ref,sc,r,smis,nstr,del,smax,smin,vmax,xc,check,flag,kis,nfld,n *g,nfg,io intent(out) vs c******************************************** call getpr(a,cov,d,dbar,e,er,f,s,v,w,mu,mn,n,nu,del,vi,va,cns,dim, *b1,b2,ro,valw,xm,pow,1,ier,io) c Error standard deviations write(io,'(/7x,''Specific Error Standard Deviations''/)') write(io,'(2x,a20,2x,f15.4)') (cand(k),dsqrt(er(k)),k=1,mu) c Eigenvalues of (2/n)*I**-1 scale=2./float(n) call eigen(a,mn,mn,f,e,s,jev) c Trap negaive f(1) if(f(1)<=0.d0) then write(io,'(/7x,''ERROR - smallest eigenvalue of max likelihood es *timates '',f12.2,'' < 0'')') f(1) write(6,'(/7x,''ERROR - smallest eigenvalue of max likelihood est *imates '',f12.2,'' < 0'')') f(1) stop endif if(f(1)<=0.d0.or.f(mn)<=0.d0) then write(io,*)'Error - negative eigenvalues from max likelihood' write(6,*)'Error - negative eigenvalues from max likelihood' stop endif semax=sqrt(scale/f(mn)) semin=sqrt(scale/f(1)) write(io,'(2x,65(''=''))') write(io,'(/7x,''Max Standard Errors of Estimates from the Eigenva *lues'',/,'' of the Covariance Matrix of the Maximum Likelihood Est *imates = '',e12.3)')semax write(io,'(/7x,''Min Standard Errors of Estimates from the Eigenva *lues'',/,'' of the Covariance Matrix of the Maximum Likelihood Est *imates = '',e12.3)')semin write(io,'(/7x,''Ratio Max/Min = '',f12.3)')semax/semin if(dim==2) then call w2(mu,v,cr,cns,r,ro,b1,b2,xm,io) elseif(dim==1) then call w1(m,v,cr,s,d,dbar,w,er,vi,vs,wv,cns,b1,ier,io) endif if(cns==1.and.dim==2) then c Constrain the mean xm(1)=xc(1) xm(2)=xc(2) endif c------- c Write errors-in-variables OLS fit parameters write(io,'(2x,65(''=''))') write(io,'(7x,''Estimates of Location & Scaling Parameters''/)') if(dim==2) then write(io,'(2x,''Ideal Point Standard Deviations = '',2(f14.2,1x)) *') b1,b2 if(vi>0) then write(io,'(/2x,''Candidate Characteristic Value Weight = '',f14. *2/)')valw endif pi=3.141592 if(cns==0.or.cns==1) then write(io,'(2x,''Rotation = '',f7.3)')180.*ro/pi elseif(cns==2) then write(io,'(2x,''Rotation = '',f7.3)')r else write(io,'(/7x,''ERROR - constraint is wrong'')') write(io,'(/7x,''cnl = '',i1)')cns stop endif c Candidate positions write(io,'(2x,65(''*''))') write(io,'(17x,''Positions of the Candidates''/)') do k=1,mu if(ref(1:1)=='y') then cr(k,2)=-1.*cr(k,2) xm(2)=-1.*xm(2) endif if(ref(2:2)=='y') then cr(k,1)=-1.*cr(k,1) xm(1)=-1.*xm(1) endif write(io,'(4x,a20,4x,f12.3,2x,f12.3)')cand(k),cr(k,1),cr(k,2) c Output candidates/parties to file 4 write(4,'(a20,1x,2(f12.3,1x))')cand(k),cr(k,1),cr(k,2) enddo write(io,'(8x,''Mean Ideal Point'',5x,f12.3,2x,f12.3)')xm(1),xm(2) endif c------ if(dim==2.and.nstr>1) then call strap(a,c,cov,d,dbar,e,er,f,s,v,w,mu,mn,n,nu,del,vi,vs,cns,d *im,bmn,br,bs1,bs2,bv,xm,pow,nstr,io) write(io,'(2x,65(''*''))') write(io,'(/15x,''Bootstrapped Parameter Means''/)') write(io,'(12x,''Rotation = '',f12.4)') 180./pi*br(1) write(io,'(7x,''Mean Ideal Point Coordinates = '',2(f12.4,2x))') *bmn(1),bmn(2) write(io,'(4x,''Ideal Point Standard Deviations = '',2(f12.4,2x)) *')bs1(1),bs2(1) if(vi>0) then write(io,'(7x,''Valence Weight = '',f12.4)')bv(1) endif write(io,'(/15x,''Bootstrapped Standard Errors''/)') write(io,'(12x,''SE of Rotation = '',f12.4)') 180./pi*br(2) write(io,'(7x,''SE of Mean Ideal Point Coordinates = '',2(f12.4,2 *x))')bmn(3),bmn(4) write(io,'(4x,''SE of Ideal Point Standard Deviations = '',2(f12. *4,2x))')bs1(2),bs2(2) if(vi>0) then write(io,'(7x,''SE of Valence Weight = '',f12.4)') bv(2) endif endif if(dim==1) then c Write parameters and candidate positions for one dimension write(io,'(/2x,''Ideal Point Standard Deviation = '',f14.2)')b1 if(vi>0) then write(io,'(2x,''Candidate Characteristic Value Weight = '',f14.2 *)')wv endif write(io,'(/17x,''Positions of the Candidates''/)') do k=1,m write(io,'(4x,a20,4x,f12.3)')cand(k),cr(k,1) enddo endif c------- c Plot ideal points if(dim==2) then call ideal2(n,nu,mu,c,uv,cr,dbar,va,valw,a,e,iw,y,vd,dav,pow,ref, *sc,io) write(io,'(''Mean of the Estimated Ideal Points'',2x,2(f12.3,2x)/ *)')y(1),y(2) call flagmed(c,rflg,m,kis,n,nu,flag,ng,check,h,med,io,4) write(io,'(2x,70(''=''))') write(4,'(''Ideal Point Mean'',5x,2(f12.3,1x))')xm(1),xm(2) write(4,'(''xbar'',17x,2(f12.3,1x))')y(1),y(2) c Voting % estimates write(io,'(7x,''% of Ideal Points that are Closest to Each Candid *ate or Party''/)') do k=1,mu if(uv(k)==1) then c Vote % write(io,'(a20,1x,f5.1)')can(k),100.*sngl(a(k))/float(nu) endif enddo endif if(dim==1) then call ideal1(nu,mu,c,cr,iw,sc,io) do k=1,nu ncode=nint(x(k,2)) write(7,'(f12.3,1x,i3)')x(k,1),ncode enddo do k=1,mu write(4,'(a20,1x,f12.3)')cand(k),cr(k,1) enddo endif c------- return end c c=========================================== c subroutine strap(a,c,cov,d,dbar,e,er,f,s,v,w,m,mn,n,nu,del,vi,va,c *ns,dim,bmn,br,bs1,bs2,bv,xm,pow,nstr,io) c******************************************** c Bootstap mean, rotation, variances using data in c c bmn - means & sigmas of mean, br - rotation, bs1 & bs2 - sigmas of x1 & x2 c w - OLS estimates (b in regres & param), xm - mean ideal point c============================================ allocatable::y real c(n,m),va(m),y(:,:),bmn(4),br(2),bs1(2),bs2(2),bv(2),xm(2) real*8 a(m,m),d(mn,mn),dbar(m),cov(mn,mn),e(m),er(m),f(mn),s(mn,mn *),v(mn,mn),w(mn*mn),s1,s2,s3,s4,rbt,z1 integer cns,dim,pow,sw,vi intent(in) c,dbar,va,m,mn,n,nu,pow,nstr,del,vi,io intent(out) bmn,br,bs1,bs2,bv c******************************************** if(nstr<2) then return endif allocate(y(0:n-1,0:m-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate work space for y in strap' stop endif s1=0.d0 s2=0.d0 s3=0.d0 s4=0.d0 do k=1,2 br(k)=0. bs1(k)=0. bs2(k)=0. bv(k)=0. enddo kount=0 do kb=1,nstr c Bootstrap c=>y call boot(c,y,e,m,n,nu,io) call getpr(a,cov,d,dbar,e,er,f,s,v,w,m,mn,n,nu,del,vi,va,cns,dim, *b1,b2,ro,valw,xm,pow,kb,ier,io) kount=kount+ier c Sum parameter estimates s1=s1+dble(xm(1)) s2=s2+dble(xm(2)) s3=s3+dble(xm(1))*dble(xm(1)) s4=s4+dble(xm(2))*dble(xm(2)) br(1)=br(1)+ro br(2)=br(2)+ro*ro bs1(1)=bs1(1)+b1 bs1(2)=bs1(2)+b1*b1 bs2(1)=bs2(1)+b2 bs2(2)=bs2(2)+b2*b2 if(vi>0) then bv(1)=bv(1)+valw bv(2)=bv(2)+valw*valw endif enddo c End of bootstrap loop if(kount==0) then write(io,'(7x,/''ERROR - each bootstrap has an error''/)') stop endif rbt=dfloat(kount) c Mean & standard deviations bmn(1)=sngl(s1/rbt) bmn(2)=sngl(s2/rbt) bmn(3)=sngl(dsqrt(s3/rbt-s1/rbt*s1/rbt)) bmn(4)=sngl(dsqrt(s4/rbt-s2/rbt*s2/rbt)) br(1)=br(1)/rbt br(2)=br(2)/rbt-br(1)*br(1) br(2)=sqrt(br(2)) bs1(1)=bs1(1)/rbt bs1(2)=bs1(2)/rbt-bs1(1)*bs1(1) bs1(2)=sqrt(bs1(2)) bs2(1)=bs2(1)/rbt bs2(2)=bs2(2)/rbt-bs2(1)*bs2(1) bs2(2)=sqrt(bs2(2)) if(vi>0) then bv(1)=bv(1)/rbt bv(2)=bv(2)/rbt-bv(1)*bv(1) bv(2)=sqrt(bv(2)) endif write(io,'(7x,''Number of Boostraps = '',i3)') kount c------- deallocate(y) return end c c=========================================== c subroutine boot(c,y,e,m,n,nu,io) c******************************************** c Bootstapped x signal in y c============================================ allocatable::w real c(n,m),y(n,m),w(:) real*8 e(m) intent(in) c,m,n,nu intent(out) e,y c******************************************** allocate(w(0:nu-1),stat=ibad) if(ibad/=0) then write(io,*) 'Unable to allocate work space for w in boot' stop endif rn=float(nu) c Uniforms call random_number(w) do k=1,m e(k)=0.d0 enddo do k=1,nu nsam=max(int(w(k-1)*rn),1) do j=1,m y(k,j)=c(nsam,j) e(j)=e(j)+dble(y(k,j)) enddo enddo c Mean score differences do k=1,m-1 e(k)=e(k)/dfloat(nu) enddo c------ deallocate(w) return end